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FOREWORD 


NASA SP-273 has been reprinted with some omissions and a few addi- 
tions to the original document, which is out of print. The reprint documents 
the current NASA Lewis chemical equilibrium computer program, CEC76. 
Since the program as well as the thermodynamic data is continuously being 
revised, the listings of these items are being omitted. However, the text is 
essentially adequate and has been reproduced without change. A few new in- 
put and output options are included in the tables for convenience. 

The following items are not being reproduced: 

(1) The program listing (appendix C) 

(2) The thermodynamic data (appendix D except pp. 168 and 169) 

(3) The examples (appendix E) 

(4) Some flow diagrams (figs. 4 and 5) 

Briefly, the additional input-output options are for 

(1) Inputting the assigned energy value on a REACTANT card in units of 

kJ/kg-mole rather than cal/g-mole. See (6a) in table VI. 

(2) Printing output tables in SI units. See SIUNIT in table VIII. 

(3) Inputting a chamber temperature estimate for rocket problems. See 

TCEST in table IX. 

(4) Freezing composition at either the throat or an exit point for rocket 

problems. See NFZ in table IX. 
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PREFACE 


This report presents the latest in a number of versions of chemical equi- 
librium and applications programs developed at the NASA Lewis Research 
Center in the past 25 years. These programs have changed over the years to 
include improved calculation techniques and to take advantage of constantly 
improving computer capabilities. 

The current version, first prepared in 1967, is based on the minimiza- 
tion of free energy approach to chemical equilibrium calculations. One pur- 
pose of this repjort is to document the program for current and future users. 
Another purpose is to present in detail a number of topics of general interest 
in complex equilibrium calculations. These topics include mathematical anal- 
ysis and techniques for obtaining chemical equilibrium; formulas for obtain- 
ing thermodynamic mixture proparties and derivatives; criteria for inclusion 
of condensed phases; calculations at a triple px>int; inclusion of ionized 
species; and applications such as constant pressure or constant volume com- 
bustion, rocket performance (to assigned area ratios or pressure ratios), 
shock wave calculations, and Chapman- Jouguet detonations . 

A number of examples are given in detail to illustrate the versatility of 
the program. 
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SUMMARY 


A detailed description of the equations and computer program for computations in- 
volving chemical equilibria in complex systems is given. A free-energy minimization 
technique is used. The program permits calculations such as (1) chemical equilibrium 
for assigned thermodynamic states (T,P), (H,P), (S,P), (T,V), (U,V), or (S,V), 

(2) theoretical rocket performance for both equilibrium and frozen compositions during 
expansion, (3) incident and reflected shock properties, and (4) Chapman -Jouguet detona- 
tion properties. The program considers condensed species as well as gaseous species. 


INTRODUCTION 

The knowledge of chemical equilibrium compositions of a chemical system permits 
one to calculate theoretical thermodynamic properties for the system. These properties 
Can be applied to a wide variety of problems in chemistry and chemical engineering. 
Some applications are the design and analysis of equipment such as compressors, tur- 
bines, nozzles, engines, shock tubes, heat exchangers, and chemical processing equip- 
ment. 

Considerable numerical calculations are necessary to obtain equilibrium composi- 
tions for complex chemical systems. This has resulted in a number of digital computer 
programs to do the calculations. (See refs. 1 and 2 for a discussion of mathematical 
procedures and references to programs.) A computer program written at NASA Lewis 
Research Center in 1961-1962 (refs. 3 and 4) has had a wide acceptance. 



Because of the still current and active interest in a general equilibrium program, it 
was considered extremely desirable to rewrite, improve, and up-date the NASA program 
(CEC71, short for Chemical Equilibrium Calculations, 1971) for the following reasons: 

(1) Many users of the previous program have made additions or deletions in order to 
do special problems. These modifications have often been difficult and/or time consum- 
ing to accomplish. CEC71 was written to make it simpler to use various parts as sub- 
routines to other programs or to include modifications. 

(2) CEC71 was written to be machine independent as much as possible. Thus, those 
parts of the previous program which were difficult to convert to other machines were 
eliminated. For example, we eliminated shift routines, subroutine BYPASS, and binary 
cards for thermodynamic data. 

(3) CEC71 incorporates programming features which have become available since the 
previous program was written. These new features include the use of namelists to sim- 
plify input, block data for common storages, and the use of logical variables and IF’s. 

(4) Additional improvements have been made in the following ways: considerable 
saving of storage, improved technique for handling condensed phases, improved method 
for obtaining assigned area ratio data in rocket performance, inclusion of a routine to 
consider discontinuities at the nozzle throat, capability of considering ionized species, 
inclusion of an option to print mole fractions of trace species, inclusion of a special pro- 
vision to handle certain singularities, inclusion of an option to calculate enthalpies of 
certain reactants, and writing chemical formulas in the output in a more conventional 
way. 

(5) A new subroutine, SHOCK, has been added to calculate incident and reflected 
shock wave parameters. 

(6) Programming has been included to permit calculation of equilibrium compositions 
with volume as one of the assigned thermodynamic functions. 

The program is now capable of doing the following kinds of problems: 

(1) Obtaining equilibrium compositions for assigned thermodynamic states. The 
thermodynamic states may be specified by the assigning two thermodynamic state func- 
tions (code names used in the program are given in parenthesis): 

(a) Temperature and pressure (TP) 

(b) Enthalpy and pressure (HP) 

(c) Entropy and pressure (SP) 

(d) Temperature and volume or density (TV) 

(e) Internal energy and volume or density (UV) 

(f) Entropy and volume or density (SV) 

(2) Theoretical rocket performance 

(3) Chapman- Jouguet detonations 

(4) Shock tube parameter calculations 
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The program may be obtained by sending a magnetic tape (minimum length, 1200 ft) 
and a written request (specifying the type of computer) to the authors. 


EQUATIONS DESCRIBING CHEMICAL EQUILIBRIUM 

Chemical equilibrium is usually described by either of two equivalent formulations - 
equilibrium constants or minimization of free energy. References 1, 5, and 6 contain 
comparisons between the two formulations. In reference 5 it was shown that, if a gener- 
alized method of solution is used, the two formulations reduce to the same number of 
iteration equations. However, in reference 6 several disadvantages of the equilibrium 
constant method are discussed. Briefly, these disadvantages are more bookkeeping, 
numerical difficulties with use of components, more difficulty in testing for presence of 
some condensed species, and more difficulty in extending the generalized method for 
nonideal equations of state. For these reasons, the free-energy minimization formula- 
tion is used. 

The condition for equilibrium may be stated in terms of any of several thermody- 
namic functions such as the minimization of the Gibbs free energy or Helmholtz free 
energy or the maximization of entropy. If one wishes to use temperature and pressure 
to characterize a thermodynamic state, the Gibbs free energy is most easily minimized 
inasmuch as temperature and pressure are its natural variables. Similarly, the 
Helmholtz free energy is most easily minimized if the thermodynamic state is charac- 
terized by temperature and volume (or density). 

Reference 1 presents equations based on minimization of Gibbs free energy. Some 
of these -equations will be repeated and expanded here for convenience. In addition, a set 
of equations based on minimization of the Helmholtz free energy will also be presented. 
However, because only ideal gases and pure condensed phases are being considered, the 
general notation of reference 1 is not used. 


Units 

In the sections prior to the COMPUTER PROGRAM section (p. 60), the International 
System of Units (SI Units) is used (ref. 7). These SI Units are as follows: 
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Physical quantity 

Unit 

Symbol 

Length 

meter 

m 

Mass 

kilogram 

kg 

Time 

second 

sec 

Temperature 

kelvin 

K 

Force 

newton 

N 

Pressure 

newton per 
square meter 

N/m 2 

Work, energy 

joule 

J 


The gas constant that is consistent with these units is R = 8314.3 J/(kg-mole)(K). In 
those sections dealing with the computer program, other units are used in addition to or 
instead of SI Units. 


Equation of State 

In this report we assume all gases to be ideal and that interactions among phases 
may be neglected. The equation of state for the mixture is 


PV = nRT 


-= nRT J 

where P is pressure (N/m 3 ), V specific volume (m 3 /kg), n moles (kg-mole/kg), 

T temperature (K), and p density (kg/m 3 ). Symbols used in this report are summa- 
rized in appendix A. Common variables used in the computer program are summarized 
in appendix B. 

Equation (1) is assumed to be correct even when small amounts of condensed species 
(up to several percent by weight) are present. In this event, the condensed species are 
assumed to occupy a negligible volume and exert a negligible pressure compared to the 
gaseous species. In the variables V, n, and p, the volume and moles refer to gases 
only while the mass is for the entire mixture including condensed species. The word 
"mixture” will be used in this report to refer to mixtures of species as distinguished 
from mixtures of reactants which will be referred to as "total reactants. " 
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The molecular weight of the mixture M is then defined to be 



(2) 


where 


n = 



1=1 


(3) 


and n. is the number of kilogram-moles of species j per kilogram of mixture. An 
equivalent expression for M is 



m 


Z 


n. 


(4) 


where M. is the molecular weight of species j. As implied in equation (4) , among the 
n possible species which may be considered, gases are indexed from 1 to m and con- 
densed species from m+1 to n. This indexing system is just to simplify the subsequent 
presentation but is not used in the program where order of species is immaterial. 


Minimization of Gibbs Free Energy 

For a mixture of n species, the Gibbs free energy per kilogram of mixture (g) is 
given by 



(5) 


where the chemical potential per kilogram-mole of species j is defined to be 
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( 6 ) 


Hi = 



The condition for chemical equilibrium is the minimization of free energy. This min- 
imization is usually subject to certain constraints such as the following mass balance 
constraints: 

n 

a ij n j “bj = 0 i = l,...,i (7) 

1=1 


or 


b.-b° = 0 i=l I (7a) 

where the stoichiometric coefficients a^ are the number of kilogram-atoms of element 
i per kilogram-mole of species j, b? is the assigned number of kilogram -atoms of ele- 
ment i per kilogram of total reactants (see eq. (191)), and 


b i ^ a ij n j 1 = 1 » ' • ‘ » 1 
1=1 


(7b) 


is the number of kilogram-atoms of element i per kilogram of mixture. 
Defining a term G to be 


l 

c=g+y\(v b °) («) 

i=l 

where X. are Lagrangian multipliers, the condition for equilibrium becomes 


+ Z Vlj) -I(V b ?)'\ - 0 (9) 

1=1 / w 
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Treating the variations 5n^ and 6X i as independent gives 

l 

+ X X i a ij = ° i = • • •» n (10) 

i=l 

and also the mass balance equation (7a) . 

Based on the assumptions in the section Equation of State (p. 4), the chemical poten- 
tial may be written 




-< 


li® + RT lnP.) + RT In P 


atm 



(j = 1 ,. . m) 

( j = m + 1 , . . . , n) 


( 11 ) 


where for gases (j = 1 to m) and for condensed phases (j > m) is the chemical po- 
tential In the standard state. For a gas, the standard state is the hypothetical ideal gas 
at unit fugacity. For a pure solid or liquid, the standard state is the substance in the 
condensed phase tinder a pressure of 1 atmosphere. The numerical values of p? that 
are generally found in the literature (ref. 8, e.g.) depend partly on a term involving 
units of atmospheres. Therefore, to be consistent, pressure P & ^ m in equation (11) 
must be In units of atmospheres. 

Equations (7) and (10) permit the determination of equilibrium compositions for ther- 
modynamic states specified by an assigned temperature T q and pressure P . That is, 
in addition to equations (7) and (10), We have the pair of trivial equations 


T = T 0 

(12a) 

P = P o 

(12b) 


However, the thermodynamic state may be specified by assigning any two state functions. 
For example, the thermodynamic state corresponding to a constant pressure combustion 
is specified, instead of by equation (12), by 


h = h 0 

(13a) 

o 

ft 

n 

ft 

(13b) 
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where h is the enthalpy of the mixture and h Q is a constant equal to the enthalpy of the 
reactants (see eq. (193)). The expression for h is 


h = 



1=1 


(14) 


where 



is the standard state enthalpy for species 


j- 


For assigned entropy and pressure (such as for an isentropic compression or expan- 
sion to a specified pressure), the thermodynamic state is specified by 


s = s 0 

(15a) 

o 

P4 

II 

P4 

(15b) 


where s is the entropy of the mixture and s Q is the assigned entropy, or entropy of the 
total reactant (see eq. (207)). The expression for s is 


where 



(16) 


(17) 


and is the standard state entropy for species j. 

The equations required to obtain composition are not all linear in the composition 
variables and therefore an iteration procedure is generally required. In the iteration 
procedure described in the following section it will be convenient to treat n as an inde- 
pendent variable. 

Gibbs iteration equations . - A descent Newton-Raphson method is used to solve for 
corrections to initial estimates of compositions n., Lagrangian multipliers X., moles n, 
and (when required) temperature T. This method involves a Taylor series expansion of 
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the appropriate equations with all terms truncated containing derivatives higher than the 
first. The correction variables used are A In n^(j = 1, . . . , m), An^(j = m + 1, . . . , n), 
A In n, = -A^/RT and A In T. As pointed out in reference 1, it is no restriction to 
start each iteration with the estimate for the Lagrangian multipliers equal to zero inas- 
much as they appear linearly in equation (10) . After making dimensionless those equa- 
tions containing thermodynamic functions, the Newton-Raphson equations obtained from 
equations (10), (7), (3), (13a), and (15a) are 



m n 

£v'j 4ln v 2 “ti a v k=1 1 <20) 

j=l j=m+l 
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(23) 

Reduced Gibbs iteration equations . - For problems with assigned thermodynamic 
states (T,P), (H,P), or (S,P), various combinations of equations (18) to (23) could be 
used to obtain corrections to estimates. However, for chemical systems containing 
many species, it would be necessary to solve a large number of simultaneous equations. 
This large number of equations can be reduced quite simply to a much smaller number 
by algebraic substitution. This is accomplished by substituting the expression for 
A In nj obtained from equation (18) into equations (20) to (23). Including equation (19) 
written with signs reversed, the resulting reduced number of equations are 



( 26 ) 
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Equations (24) to (28) are given in table I in a form which will permit direct compari- 
sons with other sets of simultaneous equations presented in later sections (tables II 
to IV). For the sake of compactness, the following symbols are used in table I: 
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J 


( 29 ) 


The script symbols subscripted with j in equation (29) are dimensionless, whereas the 
other script symbols are in units of kilogram-mole per kilogram. 

A summary of the correction equations required for several types of constant pres- 
sure problems are as follows (i = 1 to l , j = m + i to n): 


Type of problem 

Equations required 

Correction variables 

Assigned temperature 
and pressure (TP) 
Assigned enthalpy and 
pressure (HP) 
Assigned entropy and 
pressure (SP) 

(24), (25), (26) 

(24), (25), (26), (27) 
(24), (25), (26). (28) 

7T. , An,, A in n 

jt., An y A In n, A In T 

TTj, An^ A In n, A In T 


After obtaining the previous correction variables, the corrections for gaseous spec- 
ies A In n.(j = 1, . . . , m) are then obtained from equation (18). In a later section, 
Convergence (p. 24), a discussion is given on the control of the size of corrections before 
they are applied to obtain improved estimates. 


Minimization of Helmholtz Free Energy 

The equations to be presented in this section are similar to those in the previous 
section. Whatever differences appear are due to the different forms of the chemical po- 
tential pj(j = 1, . . . , m). In the previous section, pressure was one of the assigned 
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thermodynamic states and consequently the Gibbs free energy was minimized. In this 
section volume (or density) is one of the assigned thermodynamic states and consequently 
the Helmholtz free energy will be minimized. 

The two free energies have the following thermodynamic relationship: 

f = g - PV (30) 

where f is the Helmholtz free energy per kilogram of mixture, and g is given by equa- 
tion (5) - that is, 


f 


n 


= E“iV pv 

j=i 


(31) 


The chemical potential can be expressed as a thermodynamic derivative in several 
ways (ref. 9). One way is given by equation (6). Another expression is 




r 



(32) 


U 


l 

F = f +y^\.(b i -b?) (33) 

i=l 

the condition for equilibrium based on the minimization of the Helmholtz free energy sub- 
ject to mass balance constraints is 

j=i 

Treating the 6n. and 6X i as independent again gives, as in the previous section, equa- 
tions (7) and (loj. Now, however, instead of equation (11), 


‘rE ViiV"i + zl( b i- b i) ax i =0 

i=l / i=l 


(34) 
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(j=l, - • m) 


l? + RT In ^1 


n.R'T 






(j=m+l, . . n) 


( 35 ) 


where R = R/101325 to take care of units as discussed for equations (11) and (17). 

Equations (7) and (10), with pj given by equation (35), permit the determination of 
equilibrium compositions for thermodynamic states specified by an assigned temperature 
T q a °d volume V Q ; that is, in addition to equation (7) and (10), we have the pair of tri- 
vial equations 


T= T 0 

V = V Q (36) 

Analogous to equation (13) for a constant pressure combustion process, we can set 
down the following conditions for constant volume combustion: 


u' = u' 
0 

(37a) 

< 

II 

< 

o 

(37b) 


where u’ is the internal energy of the mixture and u^ is a constant equal to the internal 
energy of the reactants. The expression for u’ is 



(38) 


where is the standard state internal energy for species j. 

i 

Analogous to equation (15), for assigned entropy and volume (such as for an isen- 
tropic compression or expansion to a specified volume), the thermodynamic state is spe- 
cified by 


s = s 0 


(39a) 
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(39b) 


v = v 0 

Iteration equations are derived in the next section which permit solution of the com- 
position variables for constant volume problems. 

Helmholtz Iteration equations . - Correction equations are obtained in a manner sim- 
ilar to that described in "Gibbs iteration equations" (p. 8). In this case, however, the 
expression for Is equation (35) rather than equation (11). Because n does not ap- 
pear explicitly as a variable in equation (35), A In n no longer appears as a correction 
variable. The Newton-Raphson equations obtained from equations (10), (7), (37a), and 
(39a) are 



m n 

£ a kj n. A In n. + £ a kj A nj = bj - 1^ (k=l 1) (42) 

j=l j=m+l 
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Reduced Helmholtz iteration equations . - Equations (40) to (44) may be reduced to a 
much smaller set of working correction equations by eliminating A In n. , obtained from 
equation (40), from equations (42) to (44). Including equation (41) written with signs re- 
versed the resulting reduced number of equations are 
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x A In T 




n i (s i- R) >M 

r 2 t 


i=i 


(48) 


Equations (45) to (48) are given in table n in a form which will permit direct com- 
parison with the iteration equations in table I and the derivative equations in tables III 
and IV. For the sake of compactness the following symbols (in addition to some already 
defined in eq. (29) for table I) are used in table II: = (u^jjRT, <% = u’/RT, 

*o = u o /RT ’ and ( y v)j = (^/ R - 

A summary of the correction equations required for several types of constant vol- 
ume problems are as follows (i = 1 to l , j = m + 1 to n): 


Type of problem 

Equations required 

Correction variables 

Assigned temperature 
and volume (TV) 
Assigned internal energy 
and volume (UV) 
Assigned entropy and 
volume (SV) 

(45), (46) 

(45), (46), (47) 
(45), (46), (48) 

V 6n j 

*n., A In T 

V An jt A In T 


After obtaining the previous correction variables, the corrections for gaseous 
species A In n.(j = 1, . . . , m) are then obtained from equation (40). (See Convergence 
section (p. 24) for discussion on control of size of corrections before applying them to 
obtain improved estimates.) 

Thermodynamic Derivatives from Matrix Solutions 


All thermodynamic first derivatives can be expressed in terms of any three inde- 
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pendent first derivatives. The Bridgman tables, as tabulated for example in refer- 
ence 10, express first derivatives in terms of (8V/3T)p, (3V/3P)^,, and (3h/?T)p = c . 
We will use the logarithmic form of the volume derivatives inasmuch as in that form 
they give an indication of the extent of chemical reaction occurring among the reaction 
species. From equation (14) , 



From equation (1), 



(50) 


(51) 


Derivatives with respect to temperature . - The derivatives of n. and n with re- 
spect to temperature are needed to evaluate equations (49) and (50). These may be ob- 
tained from differentiation of equations (10), (7), and (3), which gives the following: 


l 



i-1 
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= 0 (k=l, . . o 


(54) 


m 



+ 


j=l 



j=m+l 



(55) 


As in the case of the iteration correction equations previously discussed, the deriva- 
tive equations can be reduced to a much smaller number of simultaneous equations by 
eliminating (3 In n./3 In T)p, obtained from equation (52), from equations (54) and (55). 
Including equation (53) written with signs reversed, the resulting reduced number of tem- 
perature derivative equations are 


2 m n m 



i=l j=l j=m+l j=l 



l m 



i=l j=l 




RT 


(58) 
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Equations (56) to (58) are given in table III in a form which simplifies the compari- 
son of derivative equations with the iteration equations in tables I and II. The values of 
the derivatives obtained from equations (56) to (58) could be used in equation (52) to ob- 
tain derivatives for the gaseous species (3 In n^/3 In T)p(j = 1, . . . , m), and all of the 
temperature derivatives could then be used to evaluate Cp from equation (49) . How- 
ever, there is an alternate and much simpler procedure for obtaining Cp. By substi- 
tuting (3 In n.,/3 In T)p obtained from equation (52) into equation (49) and dividing by R, 
one obtains 



In equation (59) only the temperature derivatives obtained directly from solution of 
equations (56) to (58) are required. Furthermore, all the coefficients in equation (59) are 
exactly the coefficients appearing in the reduced enthalpy equation (27). 

Derivatives with respect to pressure . - The derivative (3 In n/3 In P) T can be ob- 
tained in a manner similar to that described for obtaining derivatives with respect to the 
temperature. Differentiation of equations (10), (7), and (3) gives 


l 



i=l 


20 


ORIGINAL PAGE IS 
OP POOR QUALM 



= 0 (k=l, . . 1) 


(62) 


m 



j=l 




f? In m' 
la In P J 


n 



j=m+l 



(63) 


Equations (60) to (63) can be reduced to a smaller set by eliminating (9 In m/3 lnP) T , 
obtained from equation (60), from equations (62) and (63). When equation (61) written 
with the sign reversed is included, the results are 



Equations (64) to (66) are given in table IV for comparison with tables I, II, and III. 
The only derivative obtained from solution of equations (64) to (66) which is used is 
(3 In n/0 In P)^, (see eq. (51)). 
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Other Thermodynamic Derivatives 


As stated previously, all thermodynamic first derivatives can be expressed in 
terms of the three thermodynamic first derivatives discussed in the previous sections - 
namely, Cp, (3 In V/3 In T)p, and (3 In V/3 In P)j (see Bridgman tables in ref. 10). 
Velocity of sound a is a frequently used parameter defined by 



m 


From Bridgman tables, 



( 68 ) 


This may be written as 


where 


Using the symbols 


and 



( 68 ) 



(70) 



( 71 ) 




is 
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( 72 ) 



c 


v 


equation (69) may be written as 


y 


S ” 



(73) 


Using the equation of state given in equation (1), we obtain from equation (67) the 
familiar expression for velocity of sound 


a = 



(74) 


It should be noted that yg defined by equation (71) is required in equation (74) and not the 
specific heat ratio y defined in equation (72) . 

In the Phase Transitions and Special Derivatives section (p. 28) an alternate expres- 
sion is derived for yg for the special situation of triple phases where the expressions in 
equations (68), (69), and (73) are no longer valid. 

In reference 11 numerous first derivative relations are given that are of interest in 
rocket performance calculations. One of these relations which will be used in a later 
section is 



(75) 


PROCEDURE FOR OBTAINING EQUILIBRIUM COMPOSITIONS 

In principle, obtaining equilibrium compositions by means of the Newton -Raph son 

o 

iteration procedure discussed in a previous section should offer no difficulties. However, 
there are a number of practical items which require detailed attention in order to avoid 
numerical difficulties. These items include initial estimates, tests for condensed 
phases, phase transitions and triple points, convergence, accidental singularities, 
special handling of ions, and consideration of trace species. 
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Initial Estimates 


An extremely simple procedure is used in this report to assign initial estimates for 
composition. For the first iteration of the first point, we assign n = 0. 1. This is 
equivalent to an estimate of 10 for molecular weight. Then the number of kilogram - 
moles of each gaseous species per kilogram of mixture is set equal to 0. 1/m, where m 
is the number of gaseous species being considered. The number of moles of each con- 
densed species is set equal to zero. For HP, SP, UV, and SV problems, an arbitrary 
initial estimate of T = 3800 K is used by the program unless a different estimate is in- 
cluded with other input. 

Admittedly, this simple procedure will often give poor initial estimates. However, 
for a general chemistry program, we find this technique preferable to the alternate of 
devising numerous special routines for obtaining good estimates for numerous possible 
chemical systems. Furthermore, the estimating technique is used only for the first 
point in any schedule of points (see Main Program section, p. 79). For all points after 
the first, the results of a preceding point serve as initial estimates. 

Because no attempt is made to obtain good initial estimates, the question arises 
whether convergence can be "guaranteed". This is discussed in the next section. 


Convergence 

* 

The problem of convergence was discussed in references 1 and 3. It was pointed out 
in reference 3 that the iteration equations sometimes give large corrections, which, if 
used directly, could lead to divergence. Two situations can cause large corrections. 

The first situation occurs in the early stages of the calculation and is due to poor esti- 
mates. The second may occur at later stages of the calculation when the iteration pro- 
cess sometimes attempts to make extremely large increases in moles of species that are 
present in trace amounts. An example of this second situation is given in reference 1. 

In both of these cases a control factor X is used to restrict the size of the corrections to 
In n^(j = 1, . . . , m) and n^(j = m + 1, . . . , n) as well as to In n and In T obtained 
from solution of equations in tables I and II. 

The numerical value of, X is determined on the basis of two empirical rules which 
experience has shown to be satisfactory. For T, n, and n. for those gaseous species 

O J 

for which n./n > 10 (or ln(m/n) > -18. 420681) and for which A In n. > 0, a number 
Xj is defined as 


X, = 


1 max (|a In T|, | A In n | , | A In n. J 


(j = 1,2, . . ., m) 


(76) 
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This limits the change in T and n and the increase in n., for those species whose gas 

-8 o J 

phase mole fraction exceeds 10 , to a factor e = 7.3891. 

For those gaseous species for which ln(nj/n) £ -18.420681 and A In n. > 0, a num- 
ber Xg is defined as 

-lii^ij - 9.2103404 

X 9 = min (j = 1, . . . , m) (77) 

A In m - A In n 

-8 

This prevents a gaseous species with a mole fraction less than 10" from increasing to 
more than 10"^. The control factor X to be used in equation (79) is defined in terms of 
Xj and Xg as 

X = min(l, Xj, X 2 ) (78) 

A value for X is determined for each iteration. Whenever current estimates of 
composition and/or temperature are far from their equilibrium values, X will be less 
than 1. Whenever they are close to their equilibrium values, X will equal 1. New esti- 
mates for composition and temperature are then obtained from the following correction 
equations: 

In nj* + ^ = In nj^ + X^(A In nj)^ (j = 1, . . . , m) 

+ X^(An.)^ (j = m + 1, . . . , n) 

In n^ + ^ = In n^ + X^(A In n)^ 

In T^ i+1 ^ = In T^ + X^(A In T)^ (79) 

where the superscript i represents the i^ 1 estimate. 

The iteration procedure is continued until corrections to composition satisfy the fol- 
lowing criteria: 
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m A In 


2 -, 

j=i 


11 s 0.5xl(f 5 (j = 1,. . ., m) 


| An 
m 

I", 

j=l 


i!< 0. 5x10"^ (j = m + 1, . . . , n) 


(80) 


| A In n| < 0.5X10' 5 J 

For a constant entropy problem (SP, SV, or RKT) the following convergence test on 
entropy is also required 


s o- s 


R 


< 0.5x10 


-4 


(81) 


The convergence tests in equations (80) and (81) ensure accuracy to five places in com- 
position when expressed as mole fractions. 

In literally hundreds of different kinds of problems which have been solved by the 
CEC71 program, convergence has always been obtained in less than the 35 iterations 
permitted by the program. For the first point, which starts with arbitrary initial esti- 
mates as discussed previously, a typical number of iterations is 8 to 12. For succeed- 
ing points, which use compositions of a previously calculated point for initial estimates, 
a typical number of iterations is 3 to 5. Once in a rare while a problem may be singular 
and solutions to the correction matrix are then unobtainable. Techniques for handling 
this situation are discussed in the section Singularities (p. 30). 


Tests for Condensed Phases 

For the first point in an ENPT2 schedule of points, unless INSERT cards are used 
(see section on NAMELISTS Input (p. 64)), the program considers only gaseous species 
during iteration to convergence. For each point after the first, the program uses the re- 
sults of a previous point for its initial estimate. After every convergence, the program 
automatically checks for the inclusion or elimination of condensed species. 

The test is based on the minimization of free energy. At equilibrium, equation (9) 
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is satisfied (i.e. , 6G = 0). The requirement for a condensed species j, which was not 
previously included as a possible species, to now be included is that its inclusion will 
decrease free energy (i.e. , from eq. (9)): 



Z 


i=l 


Vij < 0 


(82) 


The subscript c in equation (82) refers to a condensed species. Equation (82) is iden- 
tical to the vapor pressure test used in reference 3 when data for the gas phase as well 
as the condensed phase are available. The vapor pressure test is 


(ifl - 


\RT 

RT \ n / 

\ 'C 



<0 


g 


(83) 


In equation (83) the subscript g refers to the gas phase of the same species as the solid 

phase referred to by the subscript c. It may be seen that equations (82) and (83) are 

l 

equivalent from the following. The term ^a.^ in equation (82) is identical for the 

i=l 

gas and condensed phases of the same species. From equation (10), this term for the 
gas phase equals (pj/HT)^ which, using the definition of in equation (11), leads di- 
rectly to equation (83) . The advantage of equation (82) over (83) is that equation (83) can 
be used only when the gas phase of the species corresponding to the condensed phase to 
be tested is present, whereas equation (82) can always be used. The use of equation (82) 
eliminates the need for the extensive programming required in reference 3 to accommo- 
date A1 2 0 3 (s, 1) for which gas phase data are not available. 

At most, only one new condensed species is included after each convergence. In the 
event that several condensed species pass the test required by equation (82), only that 
species giving the largest negative change to free energy is included as a possible species 
and convergence to a new equilibrium composition is obtained. This process is repeated 
until all condensed species required by equation (82) are included. 

If, after convergence, the concentration of a condensed species is negative, the 
species is removed from the list of currently considered species and convergence to a 
new equilibrium composition is obtained. 
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Phase Transitions and Special Derivatives 


The calculation method is based on the assumption that condensed phases are pure. 
Therefore, the possibility exists of encountering phase transition between solid and liq- 
uid (melting points) or between two stable solid phases. Such transitions constitute triple 
points since three phases of the same species coexist, one gaseous and two condensed. 
Such triple points are characterized by a definite vapor pressure and temperature, inde- 
pendent of the relative proportions of each phase. This is shown by the fact that the iter- 
ation equations of table I become singular for an assigned temperature and pressure and 
the inclusion of two condensed phases of the same species. At a triple point, for a speci- 
fied system pressure, the relative amounts of the phases can be determined only if either 
the enthalpy or the entropy is assigned. 

The program can obtain equilibrium compositions containing either one or two con- 
densed phases of a species with or without the corresponding gas phase. When tempera- 
ture is assigned, no problems arise as to which one of two or more condensed phases is 
to be considered. However, when temperature is a variable, such as in HP, SP, UV, 

SV, or RKT problems, then several possibilities exist which need to be considered. For 
example, if a liquid phase is being considered by the program and the temperature at 
convergence is below the melting point, two possibilities exist. First, the solid phase 
might be substituted for the liquid phase and a new convergence obtained, or, secondly, 
both liquid and solid might be considered simultaneously and a new convergence obtained. 
Similar possibilities exist when convergence is obtained with a solid phase above the 
melting point. These possibilities and methods of treating them are discussed in detail 
in reference 3 and the discussion will not be repeated here. In brief, the following cri- 
teria are used by the program to determine whether or not to switch one condensed phase 
of a species to another or whether to consider both simultaneously: 

Liquid present at temperature T < T m : 

If T m - T > 150 K, switch solid for liquid. 

If T m - T < 150 K, include solid and liquid. 

Solid present at temperature T > T m : 

If T - T m > 150 K, switch liquid for solid. 

If T - T m < 150 K, include liquid with solid. 

Similar tests apply for two solid phases of a species. The programming details of 
switching species are given in figure 4(c). 

The unusual situation of constant temperature during an isentropic compression or 
expansion process can occur when two phases of the same species coexist. The transition 
temperature remains constant while one phase is being converted to the other . Under this 
circumstance, derivatives with respect to temperature are not defined. Thus, several of 
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the derivatives previously discussed, Cp, c v , and (3 In V/3 In T)p, cannot be obtained 
in this case. As a consequence, equation (73) cannot be used to obtain yg. 

In reference 12, the following expression was derived for the situation of constant 
temperature and entropy: 



The velocity of sound for this situation is no longer given by equation (74) but rather by 


ag rj. — ynRTyg ip (8^) 

These derivatives will be used in connection with discontinuities in a rocket throat (see 
Discontinuities section, p. 38). 


Ions 

The program is capable of calculating equilibrium properties of plasmas (mixtures 
containing ionized species) if the plasma is considered ideal. Ideal is here meant to im- 
ply that no coulombic interactions are considered. In plasma textbooks, such as refer- 
ence 13, it is pointed out that effects of coulombic interactions do need to be considered 
in plasmas (by means of the Debye-HGckel approximation, e.g.). However, very special 
programming is needed to take care of these effects (such as accounting for ionization 
potential lowering and partition function cutoff). Therefore, the results of the plasma 
calculations made with this program will be valid only for those conditions where the 
ionic density is small so that the coulombic effects are unimportant. 

In order for ions to be considered, a charge balance equation is required: 

m 

2 a jj n j = 0 ( 86 ^ 

j=l 


where a^j indicates the excess or deficiency of electrons in the ion relative to the neu- 
tral species. For example, in a mole of an ionized species, a^ = -3 for Ar 
+1 for OZ. 


+++ 


and 


29 



The section NAMELISTS Input (p. 64) discusses the instruction required by the pro- 
gram to consider ions. 

In order to prevent difficulties in matrix solutions, the program automatically re- 
moves the charge balance equation when the mole fraction of each ion being considered 
is less than 10"®. 


Singularities 

The iteration method used in this report has successfully handled numerous chem- 
ical systems under a wide variety of thermodynamic conditions. Nevertheless, special 
procedures are required to take care of a few rare situations which would otherwise 
cause the iteration method to fail . 

One such situation is a singularity which occurs when two rows of the coefficient 
matrix are identical. This happens when the ratio of the assigned elements in these two 
rows is equal to the ratio of the stoichiometric coefficients of these two elements in every 
gaseous species for which (n./n) > 10" 8 and which contains both elements. That is, 



For example, for stoichiometric hydrogen and oxygen at 300 K and P = 1000 newtons per 
square meter, the only species with n./n > 10" 8 is HjOfe) and equation (87) applies. 
Another example is stoichiometric lithium and fluorine at 500 K and P = 101 325 newtons 
per square meter . 

A procedure which usually takes care of this difficulty is to assign values in the input 
that will make the b?/b£ ratio in equation (87) differ from the left side in the 6 th or 7 th 
figure. The program does this automatically if an equivalence ratio r = 1.0 is read in 
by reassigning r = 1.000005. Both of the previous examples are taken care of by this 
procedure. 

However, even if the b?/b£ ratio in equation (87) is made to differ from the left side 
in the 6 or 7 figure, there may be a temporary singularity. This can occur when the 
iteration procedure temporarily removes some species for which (n./n) < 10" 8 leading 
temporarily to two identical coefficient rows in the correction matrix. This difficulty is 
usually remedied by an automatic restart feature in the program. Whenever a sing ula r 
matrix occurs, the program automatically reinserts the species with n./n < 10" 8 and 
reassigns them values of nj = 10 If, after restart, the iteration procedure still leads 
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to a singular matrix, a message is printed and a test is made for another possible cause 
of singularity. 

This singularity occurs if the program is considering several condensed species 
whose determinant of stoichiometric coefficients is zero. In this case, if two or more 
of the condensed species are composed of the same chemical elements, the program will 
remove all of these species except the one with the most negative contribution to the free 
energy of the mixture based on the test of equation (82). If, after restart, a singular 
matrix again occurs, a message is printed, any results up to this point are printed, and 
the program skips to the next problem. 


THERMODYNAMIC DATA 

Thermodynamic data are included with the program for 62 reactants and 421 reaction 
species (solid, liquid, and gas phases of a species are counted as separate species). 

With few exceptions, data are taken from JANAF tables (refs. 8 and 14). Data for the 
reaction species are reduced to functional form by means of the PAC program (ref. 15). 


Assigned Enthalpies 


For each species, heats of formation (and, when applicable, heats of transition) 
were combined with sensible heats to give assigned enthalpies H^. By definition, 

5 H 298. 15 + ( H T ' ^298. 15) (88 ‘ 

We have arbitrarily assumed Hj, 8 15 - (aH?)^ ^ Equation (88) then becomes 

H T = ( aH f) 298 .15 + ( H T' H 2 88 - 15 ) 189 

In general, H° / (aH°) for T/298.I5K. For reference elements, (ah”)^ ^ 

= H° = 0. For the species included with the program these reference elements are 
Al(s)^Ar(g), B(s) (beta), Be(s), Br 2 (Z), C(s) (graphite), Cl 2 (g), Cs(s), ^(g), Fe(s), 
H 2 (g), He(g) , K(s) , Li(s) , Mg(s), N 2 (g), Na(s), Ne(g), O z (g), P(s) (red, V), S(s) (rhom- 
bic), Si(s), and Xe(g). 

Assigned enthalpies for reactants are given in table VII (in cal/mole as required for 
program input) together with some other reactant data. For cryogenic liquids, assigned 
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enthalpies are given at their boiling points . These are usually obtained by subtracting 
the following quantities from the heat of formation of the gas phase at 298.15 K: sensible 
heat between 298. 15 K and the boiling point, difference in enthalpy between ideal gas and 
real gas at the boiling point, and heat of vaporization at the boiling point. 


. Least Squares Coefficients 

For each reaction species, the thermodynamic functions specific heat, enthalpy, and 
entropy as functions of temperature are given in the form of least squares coefficients as 
follows: 



a 2 T 


a 3 T 2 + a 4 T J 


a 5 T 


(90) 


^ . a. * ^ T T 2 + T 3 * !§ T 4 ♦ 
RT 1 2 3 4 5 T 


(91) 



a, In T + a«T + ^ 1 t 2 + ^1t 3 + ?It 4 + a„ 


(92) 


For gases, the PAC program (ref. 15) first calculates thermodynamic functions from 
molecular constant data given in reference 8, and then reduces these functions to coeffi- 
cient form. For condensed species, the thermodynamic functions are taken directly from 
reference 8 for reduction to coefficient form. However, for condensed species, we first 
calculate and include data at the transition points inasmuch as these data are not given in 
the reference 8 tables. For each species, two sets of coefficients are included for two 
adjacent temperature intervals, 300 to 1000 K and 1000 to 5000 K. The data have been 
constrained to be equal at 1000 K. Details on the use of cards or tapes containing the 
coefficient data are given in section THERMO Data (p. 61). The format and listing of 
coefficient data are given in appendix D. 


CHEMICAL EQUILIBRIUM FOR ASSIGNED THERMODYNAMIC STATES 


A thermodynamic state is defined by specifying two state parameters. The following 
combination of assigned parameters are permitted by the program: 
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(1) Temperature and pressure, (T,P) 

(2) Enthalpy and pressure, (H,P) 

(3) Entropy and pressure, (S,P) 

(4) Temperature and volume (or density), (T,V) 

(5) Internal energy and volume (or density), (U,V) 

(6) Entropy and volume (or density), (S,V) 

Some of the problems handled by the program (the TP, HP, SP, TV, UV , SV, and 
SHOCK problems) use just one combination of assigned states. For example, the TP 
problem consists of a schedule of one or more assigned temperatures and pressures. 

The TP problem might be used to construct Mollier diagrams. The HP problem gives 
constant pressure combustion properties by using an assigned enthalpy equal to that of 
the reactants and a schedule of one or more pressures. Similarly, the UV problem gives 
constant volume combustion properties by using an assigned internal energy equal to that 
of the reactants and a schedule of one or more volumes. The SHOCK problem uses only 
the (T, P) combination. In appendix E (p. 187), case 950 illustrates a TP problem, case 
1565 illustrates a UV problem, and case 123 illustrates an HP problem. On the other 
hand, some of the application problems discussed in succeeding sections make use of 
more than one combination of assigned thermodynamic states. For example, the rocket 
problem uses (H,P) or (T,P) and also (S,P) and the detonation problem uses (K,P) and 
(T , P) . 


ROCKET PERFORMANCE 

The program provides routines for the calculation of theoretical rocket performance 
for both equilibrium and frozen composition during expansion. Combustion and throat 
parameters are always calculated while exit conditions are assigned optionally either as 
pressure ratios, area ratios, or both. The section DESCRIPTION OF PROGRAM INPUT 
describes in detail the various options which are permitted for this problem. Examples 
of rocket performance problems in the OUTPUT section are cases 51, 122, 679, and 

5612. 


Assumptions 

The calculation of various performance parameters are based on the following as- 
sumptions: one-dimensional form of the continuity, energy, and momentum equations; 
zero velocity in the combustion chamber; complete combustion; adiabatic combustion, 
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isentropic expansion; homogeneous mixing; ideal gas law; and zero temperature and 
velocity lags between condensed and gaseous species. For equilibrium performance, 
composition is assumed to attain equilibrium instantaneously during expansion. For fro- 
zen performance, composition is assumed to remain fixed at the combustion composition 

duA iii^ expansion. 


Parameters 

Conservation equations . - Rocket performance, as well as other fluid dynamic prob- 
lems in the program, is based on the following conservation equations: 

Continuity: 


P2 A 2 u 2 = Pl A l u l (93) 

Momentum: 

P 2 + P2 U 2 = P 1 + Pl u l 04) 

Energy: 


h 


2 




(95) 


Equation (93) describes the condition of constant mass flow rate which will be given the 
symbol m - that is, 


m = pAu (96) 

Velocity of flow . - Using subscripts c and e instead of 1 and 2 in equation (95) 
and assuming the velocity in the combustion chamber to be negligible give equation (95) 
as 


u e = V 2 < h c - h e> 


(97) 


When h is in units of joules per kilogram, u is in units of meters per second. 

Force - - From the momentum principle of fluid mechanics, the external force on a 
body in a steadily flowing fluid is due to the change of momentum of the fluid and to in- 
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crease in pressure forces acting on the body. For rocket applications, this is expressed 
as 


mu 

■ *‘ p e- p a» A e 
g c 


(98) 


The conversion factor g has been introduced to allow for various units. For some 
commonly used systems of units, such as the cgs system or International System (ref. 7), 
g c = 1. However, in the English Technical System, commonly used by engineers, 


g. 


= 32.1740 (lb mass) / JL 


(lb force) \ sec 2j 


Specific impulse . - Specific impulse is defined as force per unit mass flow rate . 
From equation (98) 


F u e, < P e- P a> A e 
rh g c m 


(99) 


In rocket literature, specific impulse is usually expressed in English engineering units 
of pounds force per (pound mass per second). However, for those systems of units pre- 
viously mentioned for which g^ = 1,1 is both dimensionally and numerically equal to 
velocity. 

In this report, when the exit pressure is equal to the ambient pressure, specific im- 
pulse will be given the symbol I g p - that is, 


I 


sp 



( 100 ) 


When the ambient pressure is assumed to be zero (vacuum conditions), specific impulse 
will be given the symbol I vac - that is, 


*vac - X sp + 


PA 
e e 

m 


( 101 ) 


Mach number . - Mach number is defined as the ratio of velocity of flow to velocity of 
sound: 
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( 102 ) 


uf = 


u 


a 


Velocity of flow is given by equation (97). Velocity of sound is given by equation (74) 

(or eq. (85)). 

Characteristic velocity . - This parameter is given the symbol c* and is defined as 


c* 



(103) 


Area per unit mass flow rate . - From equation (96) 


A = J_ 

m pu 


(104) 


The terms p and u are obtained from equations (1) and (97), respectively. 

Coefficient of thrust . - This parameter is defined in terms of previously defined 
parameters 


C 


F = 


u 

c* 


(105) 


Area ratio . - The ratio of area at any exit nozzle station to area at the throat is ob- 
tained from the values of equation (104) at these two points 


A e (A/m) e 
A t (A/m) t 


(106) 


Procedure for Obtaining Equilibrium Rocket Performance 

Rocket performance is obtained by first determining combustion properties in the 
rocket chamber and then determining exhaust properties at various stations in the nozzle 
exit including the throat. 

Combustion conditions . - Combustion temperature and equilibrium compositions are 
obtained by the program for an assigned chamber pressure and the reactant enthalpy (in 
the program HP = . TRUE . ) . From the combustion compositions and temperature the 
combustion entropy s„ and other combustion properties are determined. 
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Exit conditions . - Exit conditions may be desired for assigned pressure ratios or 
area ratios. Throat conditions, however, are always determined by the program after 
combustion conditions are completed and before any assigned pressure ratios or other 
assigned area ratios are considered. Isentropic expansion is assumed. 

For an assigned pressure ratio, equilibrium compositions and exit temperature are 
determined for the pressure P corresponding to the assigned pressure ratio and for the 
combustion entropy s c (in the program, SP = .TRUE.). For throat and other assigned 
area ratios, iteration procedures are used to determine the correct pressure ratios. 
These procedures are described in the following sections. 

After equilibrium compositions and temperature are obtained for an assigned pres- 
sure ratio or area ratio, all the rocket parameters given in the previous section may be 
determined. 

Throat conditions. - Throat conditions may be determined by locating the pressure 
ratio for which area ratio is a minimum or, equivalently, for which velocity of flow is 
equal to the velocity of sound. The second procedure is used in this report. This pres- 
sure ratio is determined by iteration. 

The initial estimate for the pressure ratio at the throat is obtained from 



^S +1 


k V^S _1) 


(107) 


The program uses the value of yg from the combustion point. 

Equilibrium properties for the P^ from equation (107) and for s c are obtained as 
for any exit point. From these properties, u^ and a? are calculated and the following 
test for convergence is made: 


u 


< 0.4x10 


-4 


(108) 


This criterion is equivalent to ensuring that, at the throat, the Mach number is within the 
value of 1±0. 2xl0“ 4 . 

If the convergence test is not met, an improved estimate of the throat pressure ratio 
is obtained from the following iteration formula: 
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(109) 
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A maximum of four iterations on throat pressure ratio is permitted. Usually, no more 
than two are required. 

Discontinuities . - A special procedure for obtaining throat conditions when velocity 
of sound is discontinuous is given in reference 12. The three steps in this procedure are 
summarized here for completeness. Steps 1 and 2 are used for all throat calculations 
while step 3 is used only for the special case of discontinuities at the throat. Case 5612 
in appendix E illustrates this particular kind of problem. 

Step 1. Obtain an initial estimate of throat pressure ratio by means of equation (107). 
For the first chamber pressure in case 5612 this gives P c /P t = 1.7469. The tem- 
perature for this estimate is 2820 K (the melting point of BeO) with liquid and solid BeO 
both present. 

Step 2. Obtain a second estimate of throat pressure ratio by means of equation 
(109). For case 5612, this gives P c /Pj. = 1.6161. The temperature for this estimate 
is 2846 K which is above the melting point with no BeO(s) present. The sequence of steps 
1 and 2 (namely, liquid and solid at one estimate of throat pressure ratio and liquid only 
at the next estimate) indicates that discontinuities exist at the throat. In this event pro- 
ceed with step 3. Otherwise, continue with the usual procedure of using equation (109) 
until convergence. 

Step 3. Estimate the pressure ratio at the melting point where the solid phase just 
begins to appear by means of 


In P. = In P + | - ? 1 -— ) (In T - In T) 
t Va In T/g m 


( 110 ) 


where the derivative is given by equation (75). The values for the right side of equation 
(110) are for the point in step 2. 

In case 5612, equation (110) gave P c /P t = 1. 7282 for which the mole fraction of 
BeO(s) equaled 0,00005. The correct value for which the mole fraction of BeO(s) is zero 
is p c / p t = 1.7277 so that the error in P c /P t using equation (110) is 0.0005 (0.03 per- 
cent). In the example in reference 12, the error is even less. 

Assigned subsonic or supersonic area ratios . - In previous programs (ref. 4), con- 
ditions for assigned area ratios were obtained by interpolation of data generated for a 
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somewhat arbitrary schedule of assigned pressure ratios. In the present program, an 
iteration scheme is used to obtain the pressure ratios corresponding to the assigned area 
ratios . 

We generally use empirical formulas that we obtained from fitting data to give initial 
estimates for the pressure ratios. However, for supersonic area ratios greater than 2, 
if the previous point were also for a supersonic area ratio greater than 2, an analytic ex- 
pression is used to obtain the initial estimate. The same analytic expression is also 
used in all cases to obtain improved estimates for the pressure ratios. 

Empirical formulas for initial estimates of P c /P e - ~ Initial estimates of pressure 

ratios corresponding to subsonic area ratios are obtained from the following empirical 
formulas: 


In — £» 


P. 


In — 


e A / A Y A e 

£+ 10.587|ln— | + 9. 454 In -£ 
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— >1.09 


( 111 ) 


and 


0.9 In 


In — £- * 


e A P / A eV A e 

— + 10. 587 (in — I + 9.454 In — 


1.0001 < — <1.09 


( 112 ) 


When an assigned supersonic area ratio requires an initial estimate of pressure ratio 
to be obtained from an empirical formula, the following formulas are used: 


P P / /A 
In — = In — + 1 / 3. 294 In — 

P e p t r \ A t 


+ 1.535 In — 


1.0001 <-£■< 2 


(113) 


and 


In — = y c + 1 . 4 In 
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(114) 
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In equation (114), the value of yg is that determined for throat conditions. The 
equilibrium properties which are obtained for the initial and each subsequently improved 
estimate of P c /P e are used in the following two equations to obtain the improved esti- 
mates of P c /P e - 

Analytic expression for improved estimates of P /P . - From table I of refer- 

' C 0 

ence 11, the following derivative may be obtained: 



1 


1 

nRT 

_ r S 

2(h c - h) 



(115) 


This derivative is used in the following correction formula to obtain an improved es- 
timate for P /P : 
c e 



where the subscript k refers to the k 411 estimate. 
The iteration procedure is continued until 




<0.00004 


(116) 


(117) 


with a maximum of 10 iterations permitted. Generally, convergence is reached within 
two to four iterations. 


Procedure for Obtaining Frozen Rocket Performance 

The procedure for obtaining rocket performance assuming composition is frozen 
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(infinitely slow reaction rates) during expansion is simpler than for that assuming equi- 
librium composition. This is due to the fact that equilibrium compositions need to be 
determined only for combustion conditions. After obtaining combustion conditions in the 
identical way described for equilibrium rocket performance, the remainder of the proce- 
dure is as follows. 

Exit conditions . - Improved estimates of the exit temperature corresponding to some 
assigned p c / p e are obtained by means of the following iteration formulas: 

( lnT 4 + r( lnT 4 + ( alnT 4 (118) 


where 



(119) 


and where k refers to the k^ estimate. The initial estimate of an exit temperature is 
the value of temperature for the preceding point. The iteration procedure is continued 
until 


| A In T | < 0.5X10" 4 (120) 

The maximum number of iterations permitted by the program is eight although the con- 
vergence criterion of equation (120) is generally reached in two to four iterations. 

Phases are also considered to be frozen. Therefore, the program will calculate 
frozen rocket performance for assigned schedules of pressure and/or area ratios only 
until an exit temperature is reached that is 50 K below the transition temperature of any 
condensed species present at the combustion point. If a calculated exit temperature is 
more than 50 K below the transition temperature, this point and all subsequent points in 
the schedule are ignored by the program and do not appear in the output. In addition, the 
following message is printed: CALCULATIONS WERE STOPPED BECAUSE NEXT POINT 
IS MORE THAN 50 DEG BELOW TEMP RANGE OF A CONDENSED SPECIES. 

After an exit temperature has been determined, all the rocket performance param- 
eters, given in a previous section, can be determined. 

Throat conditions . - Calculations for frozen throat conditions are similar to those 
discussed for equilibrium conditions. That is, equation (107) is used to get initial esti- 
mates for P /P t , equation (109) is used to get improved estimates for P c /P t , and equa- 
tion (108) is used as the convergence criterion. With composition and phases frozen, 
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there is no possibility of discontinuities at the throat as in the case of equilibrium com- 
positions . 

Thermodynamic derivatives for frozen composition . - The thermodynamic deriva- 
tives discussed in previous sections were based on the assumption that, in any thermo- 
dynamic process from one condition to another, composition reaches its equilibrium val- 
ues instantaneously. If, on the other hand, reaction times are assumed to be infinitely 
slow, composition remains fixed (frozen). In this event, expressions for derivatives be- 
come simpler (see refs. 1 and 6 for further discussion). Some derivatives, based on 
frozen composition, are as follows: 

From equations (49), (50), and (51), respectively, 



( 121 ) 


( 122 ) 


(123) 


From equation (70), 


c 


v 



- nR 


(124) 


From equations (73) and (123), it is clear that for frozen composition yg = y. 


INCIDENT AND REFLECTED SHOCKS 

The solution to the conservation equations which describe conditions for incident and 
reflected shocks is most conveniently obtained in terms of assigned temperature before 
and after shock. Therefore, theoretical values for shock parameters such as velocities 
and gas compositions are generally found tabulated as functions of assigned temperatures 
(or temperature ratios). However, in shock tube experiments, shock velocities are gen- 
erally known, rather than temperature ratios, thus requiring interpolation for purposes 
of comparison. It is therefore useful to have a calculating scheme which calculates shock 
properties in terms of assigned velocities, and that is the scheme used in this report. 
Case 52 is an example of a typical shock problem. 
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No consistent nomenclature has been found in the literature for shock parameters. 
For example, in reference 16, u is the velocity of the gases relative to the velocity of 
the shock front w, and v is the actual velocity of the gases in the shock tube. In ref- 
erence 17, however, the converse is given. We will use essentially the nomenclature 
of reference 16. 

Figure II. 6 in reference 16 shows the usual incident and reflected conditions in both 
laboratory- and shock -fixed coordinates. Velocities may be expressed as 

u = v - w = v + w or u = w - v (125) 

For incident shock waves, equation (125) gives 

Uj = Wg - Vj = Wg (126) 


and 


«2 = ^S _7 2 (127) 

In equation (126) it is assumed that the test gas is at rest (vj = 0). For reflected waves, 
equation (125) gives 


uj = Vg + Wr (128) 

and 

^5 = 7 5 + ^R = ^R (129) 

In equation (129) it is assumed that the shocked gases at the end of the shock tube are 
brought to rest (vg = 0). The asterisk is used with u£ in equation (128) to differentiate 
it from Ug in equation (127). The gas properties are the same for condition 2 but the 
relative velocities are different. All quantities that appear in equations (126) to (129) are 
positive and henceforth the arrows will be dropped. 


Incident Shocks 

The conservation equations describing incident conditions are exactly those given by 
equations (93) to (95) with area assumed to be constant (Aj = Ag). For iterative purposes, 
it is convenient to reduce these equations to a form similar to that given in reference 18 
as follows: 
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(130) 



Pl u l\/Pl 

. Pl /V 2 



(131) 


It will also be convenient to use the symbols P* and h* for the right sides of equations 
(130) and (131), respectively. These equations then become 


P* 



(132) 


h * - h 2 = 0 


(133) 


In equation (133), hg is defined by means of equation (14). 

Iterati on equations . - Applying the Newton -Raphson method to equations (132) and 

(133) divided by R and using for the independent variables the logarithm of temperature 
ratio and pressure ratio across the shock give 


3 



A 





(134) 


(135) 


where 
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and the subscript k stands for the iteration. 

The partial derivatives and right sides in equations (134) and (135) are given as 
follows: 


c 




d In V \ 
v 9 ln P /T , 2 



(137) 
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(140) 
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where, from equations (1) and (2), 


£l 

P 2 



(141) 


(142) 


(143) 


Corrections and convergence . - An arbitrary control factor X is applied to the cor- 
rections obtained from solution of equations (134) and (135) before using the corrections 
in equation (136) to obtain improved estimates. This control factor permits a maximum 
correction of 1 . 5 of the previous estimates of Pg/Pj and T 2^ T 1 • 11118 is 016 same as 
permitting a maximum absolute correction of 0. 4054652 on In Pg/Pj and In Tg/Tj. 

The control factor X is obtained by means of the following equations: 




(144) 


(145) 


X = min(X 3 ,X 4 ,l) ( 1 46 ) 

Improved estimates are then obtained by using equation (146) with equation (136) to give 
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y (147) 



The iteration process is continued until corrections obtained from equations (134) 
and (135) meet the following criteria: 




/P o\ 


A In — 

< 0.00005 

\ P l/ 


4 

< 0.00005 

w 



A maximum number of eight iterations is permitted by the program. Convergence is 

generally reached in three to five iterations. 

Initial estimates of Tg/Tj and Pg/Pj. - Formulas for temperature ratio and pres- 
sure ratio across the incident shock, assuming constant y, may be found in texts such as 
reference 16. These for mulas, slightly rearranged, are 

( 149 ) 

Pj (yj + 1) 



(150) 


If y is constant over the temperature range of Tj to 
(150) give exact theoretical values for Pg/Pj and T 2^ >r r 


Tg, then equations (149) and 
However, if y varies, then 
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equations (149) and (150) give first estimates which vary from poor to excellent. If com- 
position is assumed frozen across the incident shock, then equations (149) and (150) gen- 
erally give excellent first estimates. Improved estimates are then obtained as described 
in the previous sections. However, if equilibrium composition is assumed across the 
shock then equation (150) generally gives estimates of temperature ratio that are too 
high. An empirical correction factor is therefore used to obtain better first estimates. 

If equilibrium composition is to be assumed and equation (150) gives a value of Tg/T- 
which makes Tg > 2000 K, then the initial estimate of Tg/Tj is taken to be 

T T 

— = 0. 7 -2 (from eq. (150)) + (151) 

T 1 T 1 Tj 

If T 2 from equation (150) is greater than 2000 K, equation (151) lowers this estimate by 
three-tenths of the difference of Tg from 2000 K. 


Reflected Shocks 

The reflected shock conservation equations could be written in a manner similar to 
equations (93) to (95); however, it is more convenient to use the relations in equations 
(128) and (129) which give 


P5 W R = P2 (v 2 + W R> 

(152) 

P 5 + P5 W R = P 2 + P2( v 2 + w r) 2 

(153) 

h 5*j’4 = h 2 + j(' , 2 + w R) 2 

(154) 


For iterative purposes, it is convenient to reduce equations (152) to (154) to a form 
similar to equations (130) and (131) as follows: 
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( 156 ) 



For convenience, the right sides of equations (155) and (156) will be given the symbols 
P’ and h'. The reflected shock parameters may now be solved by the simultaneous so 

lution of 



(157) 


h’ 



(158) 


Iteration equations . - Applying the Newton-Raphson method to equation (157) and to 
equation (158) divided by R and using for independent variables the logarithm of tem- 
perature ratio and pressure ratio across the shock give 
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where 




(159) 


(160) 


(161) 
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The partial derivatives and right sides for equations (159) and (160) are given as 
follows: 




(163) 


(164) 



(165) 


(166) 
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Corrections and convergence . - The same control factors, correction equations, and 
tests for convergence discussed in the previous section for In P 2 / P l and ln T 2 ,/T l 
(eqs. (144) to (148)) are applicable for ln P g /P 2 and ln T 5^ T 2' 

Initial estimates of Tg/T 2 and Pg/Pg. - A value of T g A 2 = 2 is generally a 

satisfactory initial estimate . An estimate for Pg/P 2 in terms of Tg/Tg may be ob- 
tained from inverting equation (2.36) in reference 16. For T g /T 2 = 2 this gives 


( V si ( r S+ 1 \( P s\ 2-0 

y A^yy 


Only one solution of equation (169) is positive. 
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CHAPMAN-JOUGUET DETONATIONS 


The method used for obtaining Chapman -Jouguet detonation parameters is described 
in reference 18. There are three steps in the calculation procedure. The first step con- 
sists in obtaining an initial estimate of the detonation pressure and temperature. The 
second step is an improved estimate of these parameters by means of a recursion for- 
mula. The third step is obtaining the correct values by means of a Newton-Raphson iter- 
ation procedure. The derivation of the required equations is given in reference 18 and 
they are summarized here for convenience in slightly modified form. 

The same conservation equations (93) to (95) for continuity, momentum, and energy 
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that apply for shock also apply here with the additional constraint that u 2 = a 2> For 
iterative purposes, as was true for the shock equations, it is convenient to reduce the 
three conservation equations to the following two: 




(170) 



(171) 


The pressure ratio in equation (170) is the reciprocal of the pressure ratio in equa- 
tion (130) because, as pointed out in reference 18, the Newton-Raphson iteration en- 
counters fewer problems in the form of equation (170). 

For convenience in writing the iteration equations the symbols P" and h” will be 
used to represent the right sides of equations (170) and (171), respectively. These equa- 
tions then become 


P" 



(172) 


h” 



(173) 


Iteration Equations 


Applying the Newton-Raphson method to equation (172) and (173) divided by R and 
using for independent variables the logarithm of temperature ratio and pressure ratio 
across the detonation gives 


a 



A 





(174) 


52 



(175) 



where A In P 2 //P l and A ln T 2 //T l are defined by ec l uation ( 136 )- 

The partial derivatives appearing in equations (174) and (175) can be evaluated if y g 

is taken to be independent of temperature and pressure. This is a reasonable assumption 
for moderate ranges of temperatures and pressures. To within the accuracy of this as- 
sumption, the partial derivatives and right sides of equations (174) and (175) are 
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( 180 ) 


(h 2 - h") 
R 


( h 2 - h l) 
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y S t 2 T 2 
2M 2 




(181) 


Corrections and Convergence 

The discussion of convergence controls and corrections for shock calculations ap- 
plies equally well for Chapman -Jouguet detonations. Equations (144) to (146) give the 
control factor X; equation (147) gives the new estimates to and Tg/Tp artd 

equation (148) is the test for convergence. 


Initial Estimates of and p ^ p l 

A good initial estimate for pressure ratio is not as important as a good estimate for 
temperature ratio. For a number of chemical systems that were investigated, an initial 

estimate of fPg/PjJ =15 has been found to be satisfactory. An initial estimate for 

x / 0 

temperature ratio is found by calculating the flame temperature T 2 (HP problem) cor- 
responding to the following enthalpy: 



(182) 


I 

Improved initial estimates to the assumed value of (p^/P.) and the estimated value 

/ \ ' *'o 

of (Tj/Tj corresponding to h 2 in equation (182) can be obtained by use of the follow- 


ing recursion formulas: 
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(183) 


(184) 


where 



(185) 


(186) 


The quantities M 2 , y g 2 and (c p ) in equations (183) to (185) are the equilibrium val- 
ues for ( > 2 //p i) (^TgAij . Repeating the use of equations (183) to (186) three 

times in the program generally provides excellent initial estimates for the Newton- 
Raphson iteration. 


INPUT CALCULATIONS 

A number of options are provided in the program for specifying input and these op- 
tions are discussed in detail in the sections describing the computer program. Part of 
the input concerns reactants. The total reactant may be composed of a number of reac- 
tants and each of these reactants may be specified as an oxidant or a fuel. If the total 
reactant contains more than one oxidant, these oxidants may be combined into a total oxi- 
dant by specifying the relative proportion of each oxidant. Similarly, if the total reac- 
tant contains more than one fuel, they may be combined into a total fuel by specifying the 
relative proportion of each fuel. The overall properties of the total reactant (such as 
elemental compositions b°, enthalpy jf Qf internal energy * Q , molecular weight M q , 
density p Q , positive and negative oxidation states V + and V", specific heat <C Qi and 
entropy can then be calculated by specifying the relative amounts of total oxidant and 
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total fuel. This method is particularly convenient if calculations are to be obtained for a 
number of oxidant -fuel ratios. 

In order to obtain assigned properties for the total reactant, each reactant j must 
be specified as either a fuel or an oxidant even though a reactant may ordinarily be 
thought of as inert (as, for examole, NJ. Letting the superscript k equal 1 for oxidant 

1.L 

and 2 for fuel, the kilogram -atoms of the i in element per kilogram of total oxidant or 
total fuel is 


b i 


00 


NREAC 

z 

1=1 


a< k >n< k > 
i] ] 


i 

k = 1,2 


(187) 


where NREAC is the total number of reactants and where nj k ^ is the number of 
kilogram-moles of the j** 1 reactant per kilogram of total oxidant (k = 1) or total fuel 
(k = 2). If the amounts of oxidants and fuels are specified in terms of weights, then n| k ^ 
is obtained by 




| = 1,. . ., NREAC 
k = 1,2 


(188) 


where wj k ^ is the weight of the reactant and is the molecular weight of the j^ 1 
reactant. If, however, amounts of oxidants and fuels are specified in terms of kilogram - 
moles Nj k ^ , then 



N 

NREAC 


(k) 

i_ 



N ( k ) M ( k ) 


1=1 


j = 1 , . . . , NREAC 
k = 1,2 


(189) 


(kl 

The Mj ' may be calculated from the atomic weights of the chemical elements Mj as 
follows: 
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k = 1,2 
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i=l 


The bp) may be combined by means of the following equation to give the kilogram- 
atoms per kilogram of total reactant: 


b? = 


+ (o/f)bp) 
1 + (o/f) 


i = 1, . . . , l 


(191) 


where o/f is the oxidant to fuel weight ratio. 

Formulas analogous to equations (187) and (191) may be used to obtain other proper- 
ties of the total oxidant, total fuel, and total reactant (such as enthalpies, internal ener- 
gies, molecular weights, and densities) . 

Enthalpy of total oxidant (k = 1) or total fuel (k = 2)jT k) , (kg -mole Ag) - 


NREAC 




(k) _ 



' k)M 


,00 


RT 


k = 1,2 


(192) 


k k) 


where is the enthalpy of the 3 th reactant, (J/kg-mole)p) 

Enthalpy of total reactant kg-mole/kg: 




h o Jt + (o/f)^ 1 ) 


RT 1 + (o/f) 

Internal energy of total oxidant (k = 1) or total fuel (k = 2>l^ k ), (kg-moleAg)^'- 


(193) 
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RT 


k = 1,2 


(194) 
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where is the internal energy of the 3 th reactant, (JAg-mole)j 


(k) 


Internal energy of total reactant % Q , kg-mole/kg: 


% _ u o _ 4tr v2) + ( 0 /fXr* 1 ) 


RT 


1 + (o/f) 


(195) 


Molecular weight of total oxidant (k = 1) or total fuel (k = 2): 


M 


(k) _ 


NREAC 

I 

j=l 


(k) 


k = 1,2 


(196) 


Molecular weight of total reactant, kg/kg-mole: 


M q = 


|(o/f + iJm^m^ 
M^ + (o/f)M^ 


(If M (1) = 0, M q = M (2) and if M (2) = 0, M Q = M (1) .) 
Density of total oxidant or fuel: 


(197) 
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(198) 


Density of total reactant: 


P o = 


[.o/f) ♦ l]p<» p » 
P (1) * (o/f)p (2) 
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In the Main Program section (p. 79) several alternate expressions are given for 
specifying the relative amounts of total oxidant to total fuel and for relating them to o/f. 
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In the event that an equivalence ratio is specified, it will be necessary to make use of^ 
oxidation states. Let vt and V” be positive and negative oxidation states of the i 
element in its commonly occurring compounds. At least one of these will be zero. Thus, 
for example, the negative oxidation state for chlorine is -1 and its positive oxidation state 
is zero. The oxidation states per kilogram of total oxidant or total fuel are 



( 200 ) 


( 201 ) 


The positive and negative oxidation states for the total reactant are then 

y + _ V +(2) + (o/f)V +(1) 

1 + (o/f) 

V' = V~^ + (o/f)V~^ 

1 + (o/f) 


Equivalence ratio is now defined as 



V" 


( 202 ) 
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(204) 


One of the options in the SHOCK problem requires reactant compositions relative to 
the total reactant. These may be obtained as follows: 


(o/f)nj k) 
1 + (o/f) 


m. -4 
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(k) 


1 + (o/f) 


for k = 1 


j = l,. . NREAC 


for k = 2 


(205) 
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where m. is kilogram-mole of j^ 1 reactant per kilogram of total reactant. 

Specific heat and entropy for the total reactant in kilogram-mole per kilogram may 
be calculated using the rm, ^C°j , and Sj of the j reactants: 


NREAC 



j=l 

NREAC 



j=l 


(206) 


(207) 


COMPUTER PROGRAM 

The CEC71 computer program was written in FORTRAN IV. It was designed to 
work on any operating system compatible with FORTRAN IV. At Lewis Research Center 
it was checked out on an IBM 7094 H/7044 direct couple system, IBSYS version 13 using 
ALTIO. 

Due to shorter word length, the IBM 360 system required some double precision 
variables not required by the IBM 7094. These are indicated on the source listing in 
appendix C. Machines such as the CDC 6600 which have 60-bit words do not require any 
double precision variables. 

The source program is available to other organization as stated in the INTRODUC- 
TION. Thermodynamic data such as the data listed in appendix D as well as input for 
some sample problems (appendix E) will be included. The thermodynamic data are up- 
dated periodically. 

As a consequence, the answers for the sample problems given in appendix E may 
change somewhat from time to time. 


DESCRIPTION OF PROGRAM INPUT 

Program input data will be discussed under four categories. Three of the catagories 
are required and one is optional. The three required categories and the code names by 
which they will be referred are as follows: 
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(a) Library of thermodynamic data for reaction products (THERMO data) 

(b) Data pertaining to reactants (REACTANTS cards) 

(c) Namelist data which include the type of problem, required schedules, and options 

(NAMELISTS input) 

The optional category of data are chemical formulas of species which are singled out 
for special purposes (OMIT and INSERT cards). 

The words THERMO, REACTANTS, NAMELISTS, OMIT, and INSERT are control 
codes used by the program to identify the input read from data cards. The THERMO, 
REACTANTS, and NAMELISTS codes are punched on separate cards, each preceding 
their appropriate data cards. The OMIT and INSERT codes are on the same card as the 
data. 

The required order of the data cards is as follows: 

(1) THERMO code card 

(2) THERMO data (these cards and the preceding THERMO code card may be omitted 

if the thermodynamic data are to be read from tape) 

(3) REACTANTS code card 

(4) REACTANTS cards 

(5) OMIT and INSERT cards (if any) 

(6) NAMELISTS code card 

(7) NAMELISTS data 

In a particular rim, there may be any number of sets of REACTANTS cards (items 3 
and 4) , each followed by any number of sets of OMIT and INSERT cards (item 5) and 
NAMELISTS input (items 6 and 7). 

The program contains 3 namelists for input: INPT2, RKTINP, and SHKINP. Name- 
list INPT2 is required for all problems. Rocket and shock problems require a second 
namelist, namely RKTINP and SHKINP, respectively. 

Table V is a schematic of the input and namelist variables required for various prob- 
lems. 

Figure 1 is a flow chart of the main program which controls the reading, processing, 
and storing of input data. 

Examples of input cards for several typical problems are given in appendix E. Each 
category of input data is discussed in detail in the following sections. 


THERMO Data 

The library of thermodynamic data for reaction products are in the functional form 
discussed in the THERMODYNAMIC DATA section. The order and format of the 
THERMO data are detailed in appendix D . A listing of the card input included with the 
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program is also given in appendix D . This listing includes revisions in thermodynamic 
data through June 1970. 

THERMO data may be read either from cards or from tape. If the data are read in 
from cards, the program will write these data on tape 4. During a computer run, the 
appropriate reaction product data consistent with each new set of REACTANTS cards will 
be automatically selected from tape 4 and stored in core. 

THERMO data may be read in from cards for each run. However, a permanent tape 
or disk containing the data may be made during any run by using the required type of con- 
trol cards preceding the operating deck. 

When adding, removing, or changing data for various species on the tape, the whole 
set of THERMO data cards must be included in the input for m ak ing a new tape. 

Data for one or more species may be omitted from consideration during any particu- 
lar run (without removing the data from the tape) by using OMIT cards (see OMIT and 
INSERT cards section, p. 66). 


REACTANTS Cards 

This set of cards is required for all problems. The first card in the set contains 
the word REACTANTS punched in card columns 1 to 9. The last card in the set is blank. 
In between the first and last cards may be any number of cards up to a maximum of 15, 
one for each reactant species being considered. The cards for each reactant must give 
the chemical formula and the relative amount of the reactant. For some problems, en- 
thalpy values are required. The format and contents of the cards are summarized in 
table VI. A list of some REACTANTS cards is given in table VII. The card format is the 
same as required for the program described in references 3 and 4 . 

Relative amounts of reactants . - The relative amounts of reactants may be specified 
in several ways. They may be specified in terms of moles, mole fraction, or mole per- 
cent (by keypunching M in card column 53) or in terms of weight, weight fraction, or 
weight percent (blank in column 53). For example, in appendix E, cases 679 and 1207 
specify reactants in terms of moles and problem 51 in terms of weight. 

For these cases, the relative amounts of the reactants are completely specified by 
the values on the REACTANTS cards. However, there are optional variables which may 
be set in namelist INPT2 that indicate relative amounts of total fuel to total oxidants. 

For this situation, each reactant must be specified as a fuel or an oxidizer by keypunch- 
ing an F or O, respectively, in column 72 of the REACTANTS card. The amounts given 
on the REACTANTS cards are relative to total fuel or total oxidant rather than total reac- 
tant. 

There are four options in INPT2 for indicating relative amounts of total fuel to total 
oxidant. They include 
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(1) Equivalence ratio (ERATIO is true) 

(2) Oxidant to fuel weight ratio (OF is true) 

(3) Fuel percent by weight (FPCT is true) 

(4) Fuel to air or fuel to oxidant weight ratio (FA is true) 

For each option, the values are given in the MIX array of INPT2 (described in 
NAMELISTS Input section). This feature is illustrated by cases 52, 122, 950, and 5612 
in appendix E. Cases 52 and 950 show ERATIO = T and the reactants identified as fuel or 
oxidant in column 72. Since these cases involve just one fuel and one oxidant, the 
amounts of each (as given in columns 46 to 52) are shown as 100. This means that the 
oxidizer is 100 percent of the total oxidizers and the fuel is 100 percent of the total fuels. 
Cases 122 and 5612 are examples which have more than one fuel. Case 122 shows that 
each fuel is 50 percent (by weight) of the total fuels and the one oxidizer is 100 percent of 

the total oxidizer. 

The purpose of the previous namelist variables is to permit using one set of reactant 
cards with any number of values (up to 15) , of the variables such as oxidant to fuel ratio 
(OF = T). Case 950, for example, specifies three values of equivalence ratio 

(ERATIO = T) . 

Reactant enthalpy . - Assigned enthalpy values for initial conditions are required for 
HP, RKT, DETN, and SHOCK problems. An assigned internal energy is required for 
the UV problem. These assigned values for the total reactant are calculated automatic- 
ally by the program from the enthalpies or internal energies of the individual reactants. 

The values for the individual reactants are either keypunched on the REACTANTS 
cards or calculated from the THERMO data. The choice varies according to the type of 
problem as follows: 

(1) RKT, UV, HP problems: Enthalpies or internal energies are taken from the 
REACTANTS cards unless zeros are punched in card columns 37 and 38. For each 
REACTANTS card with the ”00" code, an enthalpy will be calculated for the species 
from the THERMO data for the temperature given in card columns 64 to 71. See MgO(s) 

in case 51, appendix E. 

(2) SHOCK problems: Enthalpies for all of the reactants are calculated from the 
THERMO data for the temperatures in the T schedule of namelist INPT2 (see table VTH) . 
If enthalpy values are punched in card columns 64 to 71 (table VI) they will be ignored. It 
is not necessary to punch zeros in card columns 37 and 38. 

(3) DETN problems: If no T schedule is given in namelist INPT2, the option for 
calculating reactant enthalpies is the same as for RKT, UV, and HP problems. How- 
ever, if a T schedule is given in INPT2, the enthalpies will be calculated from the 
THERMO data for the temperatures in the T schedule, the same as for the SHOCK 

problem. 

When the program is calculating the individual reactant enthalpy or internal energy 
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values from the THERMO data, the following two conditions are required: 

(1) The reactant must also be one of the species in the set of THERMO data. For 
example, NHg(g) is in the set of THERMO data but NHg(f) is not. Therefore, if NHj(g) 
is used as a reactant its enthalpy could be calculated automatically, but that of NH^I ) 
could not be . 

(2) The temperature T must be in the range Tj qw / 1. 2 ^ T < ^high x 1 • 2 where 
T low to T high is 1116 tem P eratur e range of the THERMO data. 

NAMELISTS Input 

As indicated in table V, the NAMELISTS code card precedes the NAMELISTS input. 
All problems require an INPT2 input. Rocket and shock problems each require an addi- 
tional set, namely RKTINP and SHK3NP. These additional sets simply follow INPT2 
directly. 

The variables in each namelist are listed in tables V, Vm, EX, and X. Table V in- 
dicates which variables are required and which are optional for the various types of 
problems. Tables Vm, EX, and X give a brief definition of each variable. Some addi- 
tional information about some of these variables follows: 

Pressure units . - The program assumes the pressure in the P schedule to be in 
units of atmospheres unless either PSIA = T, NSQM = T, or MMHG = T. 

Relative amounts of fuel(s) and oxidizer (s) . - These quantities may be specified by 
assigning 1 to 15 values for either o/f, %F, f/a, or r. If no value is assigned for any of 
these, the program assumes the relative amounts of fuel(s) and oxidizer (s) to be those 
specified on the REACTANTS cards. (See discussion in REACTANTS Cards section. ) 

Printing mole fractions of trace species . - The program automatically prints com- 
positions of species with mole fractions >5x1 0 -6 in F-format for all problems except 
SHOCK. The TRACE option permits printing smaller mole fractions. If the variable 
TRACE is set to some positive value, mole fractions greater than or equal to this value 
will be printed. When this option is used, a special E -format for mole fraction output is 
used automatically. A TRACE value of l.E - 38 is the lowest value allowed by the pro- 
gram. (See case 1565 in appendix E.) 

For SHOCK problems, mole fractions of trace species are often desired. Thus for 
SHOCK problems, the program will set TRACE to 5.E - 9 and the E -for mat for output is 
always used. This value may be changed by using the TRACE option in INPT2 namelist 
input. See case 1207 in appendix E. 

Intermediate output. - Intermediate output will be listed for whatever point IDEBUG 
is set equal to and all following points. As an example, setting IDEBUG = 3 will result 
in intermediate output for all points except the first two. 
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rpp h p. SP, TV, UV, or SV problems . - In these problems, from 1 to 26 values of 
T, P, and V (or RHO) may be assigned. However, only one value of entropy SO may be 
assigned in INPT2 for the SP or SV problem. Only one value of enthalpy is permitted for 
the HP problem and only one value of internal energy is permitted for the UV problem. 
However, these values of enthalpy and internal energy are not assigned in INPT2 but are 
calculated by the program (eqs. (192) to (195)). In a TP problem, if 26 values of T and 
26 values of P are assigned in INPT2, properties will be calculated for the 676 possible 
P and T combinations. Similarly, up to 676 combinations can be calculated for a TV 

problem. 

DETN problem. - Calculations will be made for all combinations of initial pressure 
P and initial temperature T. Initial temperatures may be specified in INPT2 namelist 
or on the REACTANTS card. 

RKT problem. - At least one chamber pressure value P is required in INPT2, al- 
though as many as 26 chamber pressures may be assigned. A complete set of calcula- 

tions will be done for each chamber pressure. 

The RKT problem requires a second namelist for input (RKTINP) discussed in the 

next section. 

RKTINP namelist (RKT problem only) . - This namelist is required for RKT prob- 
lems. It follows the INPT2 namelist. A list of variables and definitions is given m 
table IX. All variables are optional, although usually a pressure ratio schedule (PCP) , 
an area ratio schedule (SUBAR) or (SUPAR), or some combination of these schedules 

will be assigned. 

Pressure ratio and area ratio schedules must not include values for the chamber and 
throat inasmuch as these values are assigned or calculated automatically by the program. 
If both a pressure ratio schedule and an area ratio schedule are given in RKTINP, the 
pressure ratios will be calculated first. If both schedules are omitted, only chamber and 
throat conditions will be calculated. 

The program will calculate both equilibrium and frozen performance unless RKTINP 
contains FROZ = F or EQL = F. If FROZ = F, only equilibrium performance will be cal- 
culated. If EQL = F, only frozen performance will be calculated. 

SHOCK problem . - The program requires a P and T schedule in INPT2 and either 
a U1 or MACH1 schedule in a second namelist SHKINP, which is described next. These 
values all refer to the unshocked gas and must correspond one to one with each other. 

The pressure and temperature schedules are limited to 13 values for SHOCK problems 
only. This corresponds to the 13-value limit for U1 or MACH1 schedules. 

REACTANTS cards must be only for gaseous reactants that are also included as 
reaction species in the THERMO data. This permits the program to calculate enthalpy 
and specific heat values of the reactants from the THERMO data. 

The SHKINP namelist follows the INPT2 namelist. 
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SHKINP namelist (SHOCK problem only) . - A list of variables and definitions is given 
in table X. SHKINP must include from 1 to 13 values of either U1 or MACH1 of the un- 
shocked gas. The program will calculate incident shock parameters assuming both equi- 
librium and frozen composition unless SHKINP contains either INCDEQ = F, or 
INCDFZ = F. If INCDEQ = F , only frozen composition will be used. If INCDFZ = F, 
only equilibrium composition will be used. In addition there are options for calculating 
reflected shock parameters. For each incident condition called for, reflected shock 
parameters will be calculated assuming either a frozen composition (REFLFZ = T), an 
equilibrium composition (REFLEQ = T), or both (REFLFZ = T, REFLEQ = T). 


OMIT and INSERT Cards 


As indicated in table V, OMIT and/or INSERT cards may follow the REACTANTS 
cards. Their inclusion is optional. They contain the names of particular species in the 
library of thermodynamic data for the specific purposes to be discussed. Each card con- 
tains the word OMIT (in card columns 1 to 4) or INSERT (in card columns 1 to 6) and the 
names of from one to four species starting in columns 16, 31, 46, and 61. The names 
must be exactly the same as they appear in the THERMO data (appendix D) . 

OMIT cards. - These cards list species to be omitted from the THERMO data. If 
OMIT cards are not used, the program will consider as possible species all those species 
in the THERMO data which are consistent with the chemical system being considered. 
Occasionally it may be desired to specifically omit one or more species from considera- 
tions as possible species. This may be accomplished by means of OMIT cards. 

INSERT cards. - These cards contain the names of condensed species only. They 
have been included as options for two reasons. 

The first and more important reason for including the INSERT card option is that, in 
rare instances, it is impossible to obtain convergence for assigned enthalpy problems 
(HP or RKT) without the use of an INSERT card. This occurs when, by considering gases 
only , the temperature becomes extremely low (say several degrees Kelvin) . In these 
rare cases, the use of an INSERT card containing the name of the required condensed 
species will eliminate this kind of convergence difficulty. When this difficulty occurs, 
the following message is printed by the program: "LOW TEMPERATURE IMPLIES 
CONDENSED SPECIES SHOULD HAVE BEEN INCLUDED ON AN INSERT CARD" 

The second and less important reason is that if one knows that one or several partic- 
ular condensed species will be present among the final equilibrium compositions for the 
first assigned point, then a small amount of computer time can be saved by using an 
INSERT card. Those condensed species whose chemical formulas are included on an 
INSERT card will be considered by the program during the initial iterations for the first 
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assigned point. If the INSERT card were not used, only gaseous species would be con- 
sidered during the initial iterations. However, after convergence, the program would 
automatically insert the appropriate condensed species and reconverge. For all other 
assigned points the inclusion of condensed species is handled automatically by the pro- 
gram. Therefore, it usually is immaterial whether or not INSERT cards are used for the 
purpose of saving computer time. 

DESCRIPTION OF PROGRAM OUTPUT 

The program prints four kinds of output: tables of results, input data used to calcu- 
late these tables, information concerning iteration convergence, and optional interme- 
diate output. 


Tables of Results 

The final output of the program is in the form of tables that are designed to be self- 
explanatory. While each problem has its own kind of table, all tables have many features 
in common. These features are 

(1) Heading 

(2) Case number 

(3) Reactant data 

(4) Proportion of oxidant to fuel 

(5) Thermodynamic mixture properties and derivatives (P, T, p, h, s, M, 

(8 In V/0 In P)»p, (3 In V/0 In T)p, Cp, yg, and a) 

(6) Equilibrium compositon (mole fractions) 

In addition to these, certain problems list additional information: 

(1) Rocket performance: p c / p e > A e / A t» c *’ ^F’ *vac’ and *sp 

(2) Detonation: P/P p T M/Mj, p/pj, uf, and detonation velocity 

(3) Shock: ufj, Uj, u 2 , Pg/Pj, **2^1* p p^ p 1’ V 2’ u 5’ P 5 //P 2» T 5 /T 2’ 

M 5 /M 2 , p 5 /p 2 , ^d u 5 + v 2 

(4) Assigned internal energy and volume: u 

Input Data 

Input data have been previously described. The general procedure used in this pro- 
gram is to list the input as they are read in and before they are processed by the pro- 
gram. The purpose is to show, in as clear a way as possible, what is actually on the in- 
put cards. All problems list the following input data: 
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(1) The word REACTANTS 

(2) Reactant data 

(3) INSERT and/or OMIT card data 

(4) The word NAMELISTS 

(5) All data in namelist INPT2 given in table Vm (P and RHO use same storage) 
Following the INPT2 data is the statement '’SPECIES BEING CONSIDERED IN THIS 
SYSTEM". Each species in the list is preceded by some identification such as J12/65. 
The J refers to JANAF data (ref. 8). The letter L refers to unpublished data calcu- 
lated at the Lewis Research Center. The number refers to the month and the year the 
data were published or calculated (12/65 is December 1965). 

For a rocket problem, the namelist RKTINP data given in table IX are listed. For a 
shock problem, the namelist SHKINP data given in table X are listed. 

Following the list of chemical species (or RKTINP or SHKINP data, if any) is the 
current value of o/f. This is followed by a listing of the enthalpies or internal energies 
of the total fuel and oxidant (eqs. (192) and (194)) and of the total reactant (eqs. (193) and 
(195)). Following this is a list of the kilogram-atom per kilogram of each element in the 
total fuel and oxidant (eq. (187)) and in the total reactant (eq. (191)). 


Information Concerning the Iteration Procedure 

All problems except SHOCK print information regarding iteration for each point. For 
the SHOCK problem this information was omitted inasmuch as the printing of the final 
table in its present form would not be possible without extensive additional storage. Case 
122 will be used to illustrate the iteration information. Following the data on kilogram- 
atoms per kilogram of the reactants is a line containing the letters PT and the chemical 
symbols of the elements for the problem (for case 122, C, N, H, and F). The numbers 
under PT refer to the points corresponding to the columns in the final table. The num- 
bers under the chemical symbols are values of tt^. The last column (without a heading) 
is the number of iterations required to converge to equilibrium composition. In general 
there is only one line for each point unless there has been an addition or deletion of a 
condensed species. For example, in case 122 there is one line printed for point 1. For 
point 3, however, there are two lines. The first line is the result of converging with no 
condensed phases (using the estimates of the previous point) while the second line is the 
result of convergence with the inclusion of C(s) . 

For rocket and detonation problems, more than one line may be printed for condi- t 
tions other than a change in condensed species. For a rocket problem these conditions 
are for throat and for an assigned area ratio, where a line is printed out for each esti- 
mate of pressure ratio. For example, for the first point 5 of case 122, which is for an 
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assigned area ratio, the six lines indicate that six separate convergences were required 
to make a phase change and to find the correct pressure ratio for the assigned area 
ratio. For throat, additional information is given for pressure ratio and temperature 
estimates. For a detonation problem, a line is printed for each set of temperature and 

pressure estimates. 


I ntermediate Output 

The option of printing intermediate output is provided primarily as a means of ob- 
taining additional information for debugging purposes. There is usually no point in using 
this option when the program is working well. We have used this option in the past for 

the following reasons: 

(1) To find programming errors 

(2) To study the iteration process and rate of convergence 

(3) To verify that thermodynamic data have been properly prepared 

(4) To study the test for inclusion of condensed species 

As discussed in the NAMELISTS Input section intermediate output for point N and 

all following points is obtained by setting IDEBUG = N. 

Table XI is a sample of the intermediate output for the first iteration of case 123. 

It was obtained by setting IDEBUG = 1. Following the first line, which gives the iteration 
number, is the matrix corresponding to table I. The line underneath the chemical sym- 
bols contains the solution to this matrix. The next line contains some additional informa- 
tion in regard to the matrix - that is, T, n, In n (ENNL), P, ln(P/n), and X. The next 
group of lines contains information on the individual species used in setting up the pre- 
ceding matrix and the values of A In n. obtained from the matrix solution and equa- 
tion (18). The information on the individual species includes the chemical symbol, n., 

In n., dimensionless enthalpy, dimensionless entropy (SOJ/R = S°/R), negative dimen- 
sionless standard state Gibbs free enrgy (-G0J/RT £ -p°/RT), and negative dimension- 
less Gibbs free energy (-GJ/RT £ -ji./RT). 

The two derivative matrices (tables EU and IV) and their solutions are also given for 

the first point of case 123. These derivative matrices are set up after convergence, 
which, for this point, required 10 iterations. This is followed by two lines of output 
which summarize the results for point NPT. The printed variables are labeled POINT, 
P, T, H/R, S/R, MW, CP/R, DLVPT, DLVTP, GAMMA(S) , and V. The corresponding 
FORTRAN symbols which are defined in appendix B are NPT, PPP, TTT, HSUM, SSUM, 
WM, CPR, DLVPT, DLVTP, GAMMAS and VLM, respectively. 

Table XU is part of the intermediate output showing the test for condensed phases 
after convergence. The data in table XII are for case 51 (appendix E), but with a short- 
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tened pressure ratio schedule that goes directly from P c /P = 10 to P /P = 10 000. 

In appendix E it may be noted that for P c /P g = 10, which serves as the initial estimate 
for the next point, the only condensed species present is Al 2 0 3 (s). Table XII shows the 
results of testing the condensed species after convergence of P c /P = 10 000 with just 
A1 2°3^ bein S the onI y condensed species present. The notation G0-SUM(ALJ*PII) refers 
to the quantity given in equation (82). This quantity, as previously mentioned, must be 
negative for a condensed species not included in the previous convergence to be now con- 
sidered for inclusion. Further, only that species with the most negative value will be 
included. In table XII it may be seen that only one species, Ai 2 C> 3 (s), has a positive 
value for IUSE which means it was the only condensed species in the previous conver- 
gence. The temperature for this previous convergence (not shown in table XII) is 
669.9 K. For this temperature the free energy test is not made for the following species 
because their temperature range does not include this temperature: Al(Z), AlCl„(s), 
H 2 0 (s) , H 2 0(i), Mg(£), MgCl 2 (f), MgO(i), andS(s). The following species 
have positive values of G0-SUM(AU*PII) and therefore would not be considered for in- 
clusion: Al(s), AlCLjU), AlN(s), Mg(s), and S(l). The following species have negative 
values of G0-SUM(AU*PII): C(s), MgCl 2 (s) , and MgO(s) . The largest negative value for 
all preceding species tested is shown by MAX NEG DELTA G. After all condensed 
species have been tested, only that species with the largest negative value, in this case 
MgO(s) , is included as a possible reaction species and the iteration procedure is re- 
started. After convergence the next test for condensed species (not shown in table XII), 
gave no negative values for G0-SUM(AU*PII) which completes the test for condensed 
species for this point. 


Error Messages 


This section gives a list of error messages, where they occur in the program, and 
some discussion concerning them. 

CONDENSED REACTANTS NOT PERMITTED IN DETN OR SHOCK PROBLEMS 

Main program, statement 1302. As discussed in sections DETN problem and SHOCK 
oroblem (p. 65) the program accommodates only gaseous reactants for detonation or 
shock problems. After message is printed program goes on to next problem. 

CONSERVATION EQNS WERE NOT SATISFIED IN 8 ITERATIONS 

Subroutine DETON, statement 34. Convergence of conservation equations for the 
DETN problem is usually obtained in three or four iterations. The program limits the 
number of iterations to eight although we have never run a problem that required this 
many iterations. Therefore, we have not yet seen this message printed. 
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DERIVATIVE MATRIX SINGULAR 

Subroutine EQLBRM, statement 172. This situation occurs very rarely inasmuch 
as, if singularities in the matrix solutions occur, they generally occur first in the itera- 
tion matrices and the program does not get as far as the derivative matrices. However, 
it is possible for the iteration matrix to just barely avoid being singular and for the de- 
rivative matrix to be singular. When this occurs the above message is printed, DLVPT 
is set equal to -1, DLVTP is set equal to 1, and the program continues. 

DID NOT CONVERGE FOR AREA RATIO = (Value of area ratio) 

Subroutine ROCKET, statement 841. The program permits a maximum of 10 itera- 
tions to converge to the pressure ratio corresponding to the assigned area ratio. The 
usual number of iterations required is 1 to 5. The only time the number of iterations has 
exceeded 10, in our experience, has been for an assigned area ratio very close to 1 such 
as 1 . < A g /A t < 1 . 0001 . This is due to the fact that the converged throat conditions do 
not correspond exactly to an area ratio of 1 (see eq. (108)). 

DID NOT CONVERGE FOR U1 = (Value of Uj) ANSWERS PROBABLY NOT RELIABLE, 
SOLUTION MAY NOT EXIST 

Subroutine SHCK, statement 125. This message usually occurs when the assigned 
values of Uj, Tj, and Pj do not have a solution. For example, for the system contain- 
ing mole fractions of Ar = 0.9522, H2 = 0.0430, and Og = 0.0048 and for Tj = 302 K, 

Pj = 100. 2 mm Hg, the minimum value of Uj for shock to occur is 1094 meters per sec- 
ond. If Uj is assigned any value less than 1094 meters per second in SHKINP, the pre- 
vious error message will be printed. 

ERROR IN ABOVE CARD. CONTENTS IGNORED 

Main program, statement 1024. This message indicates a keypunching error in a 
control card, a missing control card, or an extraneous card. The entire contents of this 
card are ignored by the program. 

ERROR IN ORDER OF CARDS FOR (Name of species) 

Main program, statement 22. The cards containing THERMO data for the species 
named are out of order . 

ERROR IN REACTANT CARDS 

Main program, statement 52. This is due to an error in a chemical symbol such as 
the symbol not being left -adjusted or not included in BIOCK DATA. Program skips to 
next problem. 

FROZEN DID NOT CONVERGE IN 8 ITERATIONS 

Subroutine FROZEN, statement 70. Frozen expansion points generally require from 
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one to four iterations. We have not yet encountered a problem that required eight iter- 
ations. The error message was included primarily as a precaution to avoid infinite 
cycling of the iteration loop in the event of a machine error. 

INSUFFICIENT STORAGE FOR FOLLOWING (Number of) SPECIES 

Subroutine SEARCH, statement 871. This statement shows that for the chemical 
system under consideration, the program found more possible species in THERMO data 
than can be accommodated by storages reserved for the thermodynamic data in labeled 
common /SPECES/. This excess number of species is given in the error message. 

When this situation occurs, the names of the possible species are printed and control is 
returned to the main program where the next problem is read in. 

This situation can be resolved in two ways. First, the program can be recompiled 
with dimensions increased to accommodate the excess species (see section Changing 
Number of Possible Reaction Species). Secondly, OMIT cards can be used to eliminate 
the required number of excess species. 

The program is currently dimensioned for 150 species which has been found adequate 
for all problems calculated by the authors to date. 

LOW TEMPERATURE IMPLIES SPECIES SHOULD HAVE BEEN INCLUDED ON AN 
INSERT CARD, RESTART 

Subroutine EQLBRM, statement 874. This message can occur only for an HP or UV 
problem. It occurs only when the omission of an important condensed product of reaction 
causes the program to seek a combustion temperature that is unrealistically low 
(T < 100 K). When this occurs the program skips to the next problem. 

PHASES OF A CONDENSED SPECIES ARE OUT OF ORDER 

Subroutine EQLBRM, statement 156. This statement is written when the thermody- 
namic data for three or more condensed phases of a species are not in a permitted order 
as discussed in appendix D. After the statement is written, a table of output for all com- 
pleted points is written and program goes to next problem. 

REACTANT (Number) IS NOT IN THERMO DATA 

Subroutine HCALC, statement 85. As discussed in the Subroutine HCALC section 
(p. 83) enthalpies can be calculated by the program only for those reactants that are also 
included as reaction species in the thermodynamic data. The error. message is printed 
when this condition is not met. The pregram skips to the next problem. 

REACTANT TEMPERATURE OUT OF RANGE OF THERMO DATA 

Subroutine HCALC, statement 76. Subroutine HCALC permits calculation of ther- 
modynamic data for reactants for temperatures that extend up to 20 percent beyond the 
temperature range over which the data have been fitted. However, if the assigned tem- 
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perature is more than 20 percent beyond this range, the error message is printed and the 
program skips to the next problem. 

SINGULAR MATRIX 

Subroutine EQLBRM, statement 74. Singularities rarely occur in the course of ob- 
taining a matrix solution in GAUSS. When they do occur the program restarts with new 
estimates as discussed in the Singularities section (p. 30). However, if the singularity 
occurs again even with new estimates, the above error message is printed and control is 
returned to the main program which starts the next problem. One possible way to get 
through this difficulty, discussed in the Singularities section, is to assign a slightly mod- 
ified equivalence ratio or o/f. 

TEMPERATURE IS OUT OF RANGE OF THE THERMO DATA 

Subroutine SHCK, statement 1152. This message is printed whenever a converged 
temperature for a SHOCK problem is outside the following temperature range of the 
THERMO data: TLOW/1. 5 > T > 1. 25 THIGH. The program continues but the answers 
may be very inaccurate . 

THE TEMPERATURE = (degrees K) IS OUT OF RANGE FOR POINT (Number) 

Subroutine EQLBRM, statement 306. Except for the SHOCK problem, this message 
is printed whenever the converged temperature for the indicated point is outside the tem- 
perature range read in on the card following the THERMO control card. This tempera- 
ture range, which at present is 300 to 5000 K, is the one over which the gas phase ther- 
modynamic data have been fitted. Generally the thermodynamic data can be extrapolated 
a short distance without much loss in accuracy. However, to prevent large errors due to 
extrapolation if TLOW/1 . 5 > T > 1 . 25 THIGH, then after the message has been printed, 
the program writes the output for all completed points and then skips to the next problem. 

35 ITERATIONS DID NOT SATISFY CONVERGENCE REQUIREMENTS FOR THE POINT 
(Number) 

Subroutine EQLBRM, statement 973. As indicated in the Convergence section 
(p. 24), compositions are typically obtained in 3 to 12 iterations. The number 35 was 
arbitrarily selected to indicate that if convergence has not been reached by that number 
the problem probably will not converge at all. This situation occurs rarely. All com- 
pleted points up to this point are printed and program goes to the next problem. If the 
cause of nonconvergence is not obvious from the output, it is recommended that the prob- 
lem be rerun with intermediate output. 
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Sample Problems 


Nine sample problems are given that illustrate some of the features of the program. 
Four are rocket performance problems, RKT = T, (cases 51, 122, 679, and 5612); two 
are combustion problems (case 123 is for combustion at constant pressure, HP = T, and 
case 1565 is for combustion at constant volume, UV = T); case 52 is a detonation prob- 
lem, DETN = T; case 1207 is a shock problem, SHOCK = T; and case 950 is an as- 
signed temperature and pressure problem, TP = T. 

It would not be practical to illustrate every possible variation of options permitted 
by the program. However, the sample problems were selected to illustrate many of the 
possible variations and in particular those variations which we feel might most often be 
used. Included in the features illustrated are the following: 

(1) Specifying proportions of various reactants 

(a) o/f: cases 122, 123, 1565 

(b) Equivalence ratios: cases 52, 950 

(c) Percent fuel by weight: case 5612 

(d) Complete information on reactant cards: cases 51, 679, 1207 

(e) Relative weights of reactants: cases 51, 52, 122, 123, 950, 1565, 5612 

(f) Relative moles of reactants: cases 679, 1207 

(2) Specifying enthalpies 

(a) On reactant cards: cases 51 (partly), 122, 123, 679, 950, 1565 (partly), 

5612 

(b) Calculated by program: cases 51 (partly), 52, 1207, 1565 (partly) 

(3) Pressure units 

(a) psia: cases 51, 122, 679 , 56 1 2 

(b) atm: cases 52, 123, 950 

(c) mm Hg: case 1207 

(4) INSERT: cases 51, 5612 

(5) OMIT: case 950 

(6) Composition in floating point format: case 1565 

(7) Program considers ions: case 679 

(8) Special derivatives due to two condensed phases of a species: cases 51, 5612 

(9) Special throat interpolation: case 5612 

Some additional features of the program illustrated by the various cases are the 
following: 

Case 51: This case shows several condensed species being automatically inserted 
and removed by the program. Frozen expansion is stopped at point 3 inasmuch as exit 
temperature is below melting point of 2315 K. 
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Case 122: This case shows that it is possible to assign a schedule of points which 
include a mixture of pressure ratios, subsonic area ratios, and supersonic area ratios 


MODULAR FORM OF THE PROGRAM 


In order to facilitate adding or deleting applications of the chemical equilibrium part 
of the program, the program was set up in eight modules. These modules are concerned 
with general input, additional input processing, four applications, equilibrium calcula- 
tions, and output. The general flow of these modules and associated routines are given 
in the following schematic: 



A subroutine tree diagram is given in figure 2. From these diagrams it is clear that, 
for example, the rocket application could be eliminated by omitting subroutines ROCKET 
RKTOUT, and FROZEN and by omitting the statement which calls ROCKET in the main 


program. 

Some details of the individual routines are given in 


the ROUTINES section (p. 78). 
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The following section gives: (1) the general purpose of each module, and (2) common 
variables set by the module which are either required or set by the equilibrium module. 


General Input Module 

This module is controlled by the main program. A flow diagram of the main pro- 
gram is given in figure 1 . The module sets up input required by all application modules. 
Some of these input data include the following: 

(1) THERMO data. The main program reads the thermodynamic data cards for all 
species and writes the data on tape 4. Subroutine SEARCH pulls the data for the appro- 
priate species for a particular chemical system from tape 4 and stores the data in core. 

(2) REACTANTS cards. These cards are read and processed by subroutine REACT. 

(3) INPT2 namelist data. The main program initializes the variables, reads and 
writes the data, and converts some of the data to the form required by the remainder of 
the program. 

(4) BLOCK DATA. These data are set and remain for the entire computer run. 

(5) Composition estimates. These estimates and associated variables are set for the 
first iteration for the first point only. 

The variables having to do with composition which are used directly by the equilib- 
rium module include ENN, ENNL, IQl, IUSE, JSOL, JUQ, EN(j,l), and ENLN(j) . (All 
common variables are defined in appendix B.) The remaining variables set in the general 
input module which are used by the equilibrium module include A, COEF, IONS LLMT 
RR, NS, NLM, NC, SHOCK, IDEBUG, TEMP, TRACE, and TMID. These latter varial 
bles usually remain unchanged for an entire problem set. 


Application Modules 

Each application module is called from the main program according to the type of 
problem designated in namelist INPT2. The module controls the flow of the program 
until the problem is completed. Control is then returned to the main program where the 
next code card is read. A general flow diagram is given in figure 3. The module is re- 
sponsible for the following items: 

(1) Reading any additional input for the particular problem (namelist RKTENP for the 
RKT problem and SHKINP for the SHOCK problem) . 

(2) Doing any calculation peculiar to the problem. 

(3) Calling routines in the additional input processing module (p. 75) as required. 
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(4) Setting a point number NPT and setting certain variables required for the as- 
signed thermodynamic states for that point if they have not already been set in the main 
program. These variables are set as follows for various assigned states: 






Variable 

Temperature 
and pressure 

Enthalpy and 
pressure 

Entropy and 
pressure 

Temperature 
and volume 

Internal energy 
and volume 

Entropy and 
volume 

TP 

HP 

SP 

VOL 

TT 

PP 

VLM (NPT) 
SO 

HSUBO 

True 

False 

False 

False 

T o 

P o 

False 

True 

False 

False 

Estimated T 
P o 

h c /R 

False 

False 

True 

False 

Estimated T 
P o 

Jq 

True 

False 

False 

True 

T o 

V o 

False 

True 

False 

True 

Estimated T 
V o 

False 

False 

True 

True 

Estimated T 
V o 

Jo 


Estimated T is set equal to 3800 K for the first point only . For succeeding points the 
temperature estimate is set equal to the value of some previous converged point. 

(5) Calling the equilibrium module (subroutine EQLBRM). 

(6) Printing special output for the problem and calling the output module to print 
general output. 


Equilibrium Module 


The equilibrium module calculates compositions and thermodynamic properties for a 
particular point NPT. The module is controlled by subroutine EQLBRM which is dia- 
grammed in figures 4(a), (b), and (c). 

Subroutine EQLBRM calls three subroutines: (1) CPHS to calculate thermodynamic 
functions of the individual species, (2) MATRIX to set up the matrix according to tables 
I n m , and IV as required, and (3) GAUSS to solve the set of equations. 

Common variables either calculated or set by the EQLBRM module include CPR, 
CPSUM, DLVPT, DLVTP, EN, ENLN, ENN, GAMMAS, HO, HSUM, JLIQ, JSOL, PPP, 
S, SSUM, TOTN, TTT, VLM, andWM. 


77 




Additional Input Processing Module 

The additional input processing module consists of a set of routines which are called 
for various purposes: 

(1) NEWOF is called by an application module to adjust values of b? D r and 

u q/R or h Q /R for each o/f. 1 ’ °’ 

(2) SAVE is called to save or move composition data from one point to another. The 

purpose is to use the calculated results from some previous point as initial estimates for 
the current point. 

(3) HCALC is called from three places: 

NEWOF - to calculate enthalpies for reactants where one or more REACTANTS 
cards had zeros in columns 37 and 38 

SHCK - to calculate compositions and thermodynamic properties when compo- 
sition is frozen to the reactants composition 

DETON - to calculate enthalpy, HSUBO, for reactants at the temperatures as- 
signed in the T array in namelist INPT2 

(4) CPHS is called by HCALC to calculate thermodynamic functions for an individual 
reactant using the THERMO data. 

Rn by 11118 module * hich are required by the equilibrium module include 

BO, HSUBO, and EQRAT. 


Output Module 

The output module consists of the three subroutines VARFMT, EFMT and OUT1 with 
two entrlee OUT2 and OOT3. OUT1 lisle data given on the REACTANTS cards as well 
as o/f, %F, r, and p . OUT2 lists the properties P . , T o(k/cc) hfcal/e) 

s<ca./,g>«K», ,3 In V/3 In P, T . ,3 In V/3 In T) p , lists 

equilibrium mole fractions of the reaction species 

Subroutines VARFMT and EFMT are called from OUT2 and OUT3. VARFMT ad- 
justs the number of decimal places in a variable format according to the size of the num- 
ers. EFMT sets up a special E-type format for printing density p and mole fractions. 


ROUTINES 


SSESSSS 


The program consists of a main program, 17 subroutines, 3 entries, and block data 
A description of the function of each of these is given in the following sections 

Most program variables used in these sections are in labeled common and are de- 
fined m appendix B. 
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Main Program 


A flow chart of the main program is given in figure 1. Generally , the routine per- 
forms the following functions: 

(1) Reads code cards THERMO, REACTANTS, OMIT, INSERT, and NAMELISTS 

and directs flow of program accordingly. 

(2) Stores THERMO data on tape 4. 

(3) Calls subroutine REACT to read and process REACTANTS cards. 

(4) Reads OMIT and INSERT cards and stores species names. 

(5) Initializes variables in namelist INPT2. 

(6) Reads and writes namelist INPT2. 

(7) Converts assigned densities, if any, (RHO(i) in INPT2) to specific volumes: 
VLM(i) = l/RHO(i). 

(8) Stores the number of pressures or volumes in NP. 

(9) Stores values of o/f inOXF array. If o/f values have not been inputted direc- 
tly, they are calculated as follows: 


Values In MIX (table VIII) 

Code (table VIII) 

,o/f calculation in main program 

Oxidant to fuel weight ratio, o/f 
Fuel to air weight ratio, f/a 
Percent fuel, %F 

Equivalence ratio, r 
Not specified 

OF = .TRUE. 

FA = .TRUE. 
FPCT = .TRUE. 

ERATIO = . TRUE . 

OXF(i) = o/f 
OXF(i) = l/(f/a) 

OXF(i) = (100-%F)./(%F) 

-rV"^ - V +(2) 
OXF(i) = — 

r V" (1) + V* 111 

OXF(i) . WP »> 

WP(2) 


Values of WP(1) and WP(2) are defined in appendix B. 

(10) Makes necessary adjustments to consider charge balance if IONS - .TRUE. . 

This is done by adding 1 to NLM and E to LLMT array. 

(11) Calls SEARCH to pull required THERMO data from tape and to store the data in 

core . 

(12) Sets initial estimates for compositions. These estimates are set with each 
INPT2 read. They are used only for the first point in the lists of variables in namelist 
(e.g., the first o/f and the first T and P in a TP problem) . All succeeding points 

use results from a previous point for estimates. 

For the first point the program assigns an estimate of 0. 1 for n, the total number of 
kilogram -moles per kilogram. The initial estimate of number of moles of each gaseous 
species per kilogram of mixture n. is set equal to 0. 1/m where m is the total number 
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of gaseous species. Condensed species are assigned zero moles. (See Initial Estimates, 
p. 24)). 

(13) Sets IUSE(j) positive for condensed species listed on INSERT cards (see IUSE 
array) . 

(14) Calls THERMP if TP, HP, SP, TV, UV, or SV is true. 

(15) Calls DETON if DETN is true. 

(16) Calls SHCK if SHOCK is true. 

(17) Calls ROCKET if RKT is true. 


Subroutine REACT 

The purpose of subroutine REACT is to read and process the data on the REACTANTS 
cards . The subroutine is called from the main program after a REACTANTS code card 
has been read. The data on these cards are described in the REACTANTS Cards section 
(p. 62). 

The reactants may be divided into two groups according to card column 72 on the 
REACTANTS cards. The two groups are oxidants (O in cc 72) and fuels (cc 72 /O). We 
generally keypunch F in card column 72 for fuels even though this is not necessary. The 
contents of card column 72 are read into FCX. Depending on the contents of FOX, pro- 
gram variables relating to oxidants or fuels are subscripted 1 for oxidants and 2 for fuels . 

The FORTRAN symbols for the properties read from the REACTANTS cards and 
their associated properties discussed in INPUT CALCULATIONS (p. 55) are as follows: 


Property 

FORTRAN symbol 

a (k) 

a ij 

ANUM(j,m) a 

w (k > 

J 

PECWT(j) (if no M in cc 53) 

n!» 

J 

PECWT(j) (if M in cc 53) 

(<’ 

ENTH(j) (if not UV problem and 00 not in cc 37 and 38) 

(< 

ENTH(j) (if UV problem and 00 not in cc 37 and 38) 

<•? 

DENS(j) 


a Each of the j REACTANTS cards contains from 1 to 5 stoich- 


iometric coefficients read (indicated by subscript m) into 
ANUM(j , m) and their corresponding chemical symbols read 
into NAME(j.m). In relating an ANUM(j,m) with the 

index i associated with a particular chemical element is 

determined from the chemical symbol in NAME(j,m). G 

ORIGINAL PAGE'S 
OF POOR QUALD-i 




If there are several oxidants their properties are combined by subroutine REACT 
into properties of a total oxidant using the relative proportion of each oxidant given on the 
REACTANTS cards. Similarly, if there are several fuels, their properties are com- 
bined into properties of a total fuel. The total oxidant and total fuel properties discussed 
in INPUT CALCULATIONS and their associated FORTRAN symbols are as follows. 


Property 

FORTRAN symbol 

Equation 

b [ k > 


B0P(i,k) : 

(187) 

M< k > 


RMW(j) 

(190) 

jtMt 


HPP(k) (if not UV problem and 

(192) 



00 not in cc 37 and 38) 




HPP(k) (if UV problem and 00 

(194) 



not in cc 37 and 38) 




AM(k) 

(196) 

p (k) 


RH(k) 

(198) 

v+00 


VPLS(k) 

(200) 

/ 1 

'JREAC 

VMIN(k) 

(201) 

wj k y 

2 ^ wj k) 

PECWT(j) 


/ 

H 




If any of the are zero then RH(1) = RH(2) = 0. 

These total oxidant and total fuel properties are subsequently combined into total re- 
actant properties by using the values of oxidant-fuel mixture ratios obtained from the 
main program. This is done in NEWOF, an entry in SAVE. 

Other common variables set by REACT are LLMT, NAME, ANUM, ENTH, FAZ, 

RTEMP, FOX, DENS, RMW, MOLES, NLM, NEWR, and NREAC. 

A provision is made for eliminating a second tape search when two consecutive sets 
of REACTANTS cards contain the same elements. This is done by saving the element 
symbols (LLMT(Z)) in LLMTS(Z), the kilogram -atoms per kilogram (B0P(Z,k)) m 
SB0P(Z ,k), and the number of elements (NLM) in NLS. 

Atomic weights used in equation (190) are stored in ATOM(2,i). The corre- 
sponding chemical symbols are stored in ATOM(l.i). The oxidation states of the chemi- 
cal elements V+ or V: used in equations (200) and (201) are stored in ATOM(3,i). The 

ATOM array is stored by BLOCK DATA. 
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Subroutine SEARCH 


Subroutine SEARCH reads the THERMO data which have been stored on tape 4 and 
stores the appropriate data in core . 

A check is made near the beginning of the routine to prevent THERMO data from ex- 
ceeding their storage allotments These variables are ail in labeled common SPECIES 
and are currently dimensioned for 150 species (appendix B) However, this dimension 
may be reduced to save storage (see Changing Numer of Possible Reaction Species, 

P- 91). 

SEARCH is called from the main program when the logical variable NEWR is true 
(see fig. 1). NEWR is set true in REACT to indicate a new chemical system. REACT 
also stores chemical element symbols for the current chemical system in the LLMT 
array. SEARCH stores THERMO data in core for each species whose elements are in- 
cluded in the LLMT array (unless the species name was listed on an OMIT card). 

The THERMO data are stored in common variables TLOW, TMID, THIGH, SUB, A, 

and TEMP. SEARCH writes out the names and dates of species whose data are 
stored in core. 

SEARCH initializes the IUSE array. IUSE(j) f.or gaseous species are set equal to 
zero. IUSE(j) for condensed species are set equal to negative integers. For the chemi- 
cal system under consideration, the first possible condensed species is set equal to -1 , 
the second to -2, and so on, with one exception. In the event there are two or more con- 
densed phases of the same species, each phase is given the same negative integer. Thus, 
if IUSE(j) for B 2 Og(f) is set equal to -4, for example, IUSE(j) for B 2 0 3 (s) will also be 
set equal to -4. A description of the IUSE array is given in the next section. 

The various condensed phases of a species are expected to be adjacent in the 
THERMO data as they are read from tape 4. These phases must be either in increasing 
or decreasing order according to their temperature intervals. 

NS contains the total number of species stored in core. NC contains the total number 
of condensed species (counting each condensed phase of a species as a separate species) . 

IUSE array. - Each value in the IUSE array is associated with a species. These val- 
ues of IUSE serve two purposes: 

(1) They indicate which species are to be included in the current iteration (IUSE(j) 

< 0 for excluded species and IUSE(j) > 0 for included species). 

(2) They indicate multiple phases of the same species if absolute values of IUSE(j) 
are equal. 

The IUSE(j) are initialized in subroutine SEARCH and the main program as follows: 

(1) IUSE(j) = 0 for all gaseous species. 

(2) IUSE(j) = n for all condensed species whose names have been listed on INSERT 
cards. The number n indicates the species was the n th condensed species whose 
THERMO data were read from tape 4. 
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(3) IUSE(j) = -n for all condensed species not listed on INSERT cards where n is 

dMn rZ ( iLl values of .USE(j) may be adiusted later in subroutine EQLBRM. For 
condensed species, the sign is adjusted as species are included or excluded m the cur- 

1:601 For tiie IONS option, IUSE(j) values for ionic secies are set to -10000 when the 
mole fractions of all ionic species are less than 10 . 


Subroutine HCALC 

The purpose of HCALC is to calculate thermodynamic properties to reactonts under 
certain circumstances. HCALC is called from NEWOF (see "Entry NEWOF ), SHCK, 

^ is called from NEWOF when CALCH is set true. CALCH is set true in the 

main program when zeros have been punched in card columns 37 and 38 on one or more 
REACTANTS cards. The zeros are a code indicating that the enthnlps ^ at 

ergy for UV problems) for the reactant should be calculated from the THERMO data at 
the temperature punched on the card. This temperature has been stored m the R 

array. CPHS is called to calculate the enthalpy . The value is stored in the ENTH ar y 

and printed in the final tables fortran symbols, and the 

The properties calculated in subroutine HCALC, 

conditions for which they are used are as follows: 


Property 

FORTRAN symbol 

Equation 

Conditions 


HPP(k) 

(192) 

SHOCK problem. DETN problem with T schedule. 




HP, RKT, or DETN problem if 00 in cc 37 and 38 

h /R 

HSUBO 

(193) 

SHOCK problem. DETN problem with T schedule. 

n 0 / 



HP, RKT, or DETN problem if 00 in cc 37 and 38 

*(k) T 

HPP(k) 

(194) 

UV problem If 00 in cc 37 and 38 

u’ /R 

HSUBO 

(195) 

UV problem if 00 in cc 37 and 38 

M 

AMI 

(197) 

SHOCK or DETN problem 

0 

m 

EN(j) 

(205) 

SHOCK problem 

i 

CPR1 

(206) 

SHOCK problem 

*0 

SO 

(207) 

SHOCK problem 


The quantity m* was deliberately subscripted differently from EN(j) to allow for the 

fact that the same compound may have a different index as a reac ^^ n ™ ^ ^CTANTS 
species. Thus, for example, 0 2 (g) might be the third reactant read in from REACTA 
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cards and also the tenth species read in by SEARCH. In this case m, would be stored 
in EN(10) . 


Subroutine SAVE 

This subroutine has several functions, all of which are concerned with saving some 
information from a completed calculation for subsequent use in later calculations . The 
primary purpose is to save computer time by having good initial estimates for composi- 
tions. 

These estimates for the next point, NPT, come from either the point just completed, 
ISV , or some other previous point. The flow of the routine is directed by ISV as follows: 

(1) ISV positive. Transfer compositions for point just completed for use as initial 
estimates for next point (transfer EN(j,ISV) to EN(j,NPT)). 

(2) ISV negative. Save values of ENLN(j) for gases and EN(j) for condensed in 
SLN(j), ENN in ENSAVE, ENNL in ENLSAV, IQl in IQSAVE, JSOL in JSOUS, JLIQ in 
JLIQS, and NLM in LL1. (These values are saved because they are to be used as initial 
estimates for some future point and they may be overwritten in the meantime.) Make 
ISV positive and transfer EN(j,ISV) to EN(j,NPT). 

(3) ISV zero. Use the data previously saved (as discussed in (2)) as initial estimates 
for current point. Restore IUSE codes and inclusion or exclusion of "E M as an element 
for IONS option. 

Entry NEWOF . - NEWOF combines the properties of total oxidant and total fuel cal- 
culated in subroutine REACT with an o/f value to give properties for the total reactant. 
NEWOF is called for each mixture assigned in the MIX array in INPT2 namelist. It is 
called from either THERMP , ROCKET, SHCK, or DETON. The calculated properties 
and corresponding FORTRAN symbols are as follows: 


Property 

FORTRAN symbol 

Equation 

b i 

B0(i) 

(191) 

h 0 /R 

HSUBO (if not UV problem) 

(193) 

u^/R 

HSUBO (if UV problem) 

(195) 

P 0 

RHOP 

(199) 

r 

EQRAT 

(204) 




Subroutine HCALC is called by Entry NEWOF to calculate the enthalpies for each re- 
actant that has zeros keypunched in card columns 37 and 38 in its REACTANTS card (see 
table VI) . 

Values of KPP(2), KPP(l), HSUBO, B0P(i,2), B0P(i,l), and B0(i) are printed out. 
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Subroutine EQLBRM 


EQLBRM is the control routine for the equilibrium module which calculates equilib- 
rium compositions and thermodynamic properties for a particular point. The subroutine 
flow diagram is given in figures 4(a) to (c) . Figure 4(a) is a diagram of the complete 
routine. Figure 4(b) is an expansion of the block labeled 85 in figure 4(a). It gives the 
details for obtaining and applying corrections in subroutine EQLBRM. Figure 4(c) is the 
expansion of the block between blocks labeled 160 and 143 in figure 4(a) . Figure 4(c) is 
the flow diagram for adding and removing condensed species. 

Several common variables must be set before EQLBRM is called. These variables 

and the modules setting them are summarized as follows: 


Module 

Variables 

Input 

ENN, ENNL, IQ1, IUSE(j), JSOL. j for first iteration. 
JLIQ, EN(j.NPT), ENLN(j) J first point only 

COEF(k,i,j), TEMP(j.k), IONS. LIMT(Z), RR, NS 
NLM, NC, SHOCK, IDEBUG, TRACE, and TMID. 


For TP, HP, SP, TV, UV, and SV problems, TP, HP, 
SP, VOL, and SO (if SP is true) are all set. 

Application 

TT, PP, NPT, VLM(NPT) (if VOL is true). In addi- 
tion, if not previously set in Input: VOL, TP, HP, 
SP, and SO (if SP is true). 

Additional input 

BO, EQRAT, and HSUBO 

j processing 



Common variables set by EQLBRM include TTT(NPT), PPP(NPT), SSUM(NPT), 
HSUM(NPT), CPR(NPT), GAMMAS (NPT) , VLM(NPT), WM(NPT), DLVPT(NPT) , 
DLVTP(NPT) , TOTN(NPT) , ENN, EN(j,NPT), ENLN(j), IUSE(j), JLIQ, and JSOL. 


Subroutine CPHS 

Subroutine CPHS calculates thermodynamic properties using equations (90), (91), 
and (92) for species numbering from 1 to NS for an assigned temperature TT. It uses 
either one of two sets of coefficients: COEF(2, i,j) for the temperature interval TLOW 
to TMID and COEF(l,i, j) for the interval TMID to THIGH. The index j is the the j 
species and the index i (i = 1 to 7) refers to the 1 th coefficient. 

The properties calculated and their corresponding FORTRAN symbols are as 

follows: 



Property 

FORTRAN symbol 

{ s t / r ). 

S(j) j = JSl,.. 

., NS 


K0(j) j = JSl,. 

. NS 

NS 

j=JSl 

CPSUM (one of the terms in eq. (59)) 

"j( c ?/ r ), 

CPSUM (when CPHS is called by HCALC) 


The index JS1 is always set equal to 1 in all routines calling CPHS except HCALC. In 
the latter event JS1 and NS are both set equal to the index j of a particular species. 


Subroutine MATRIX 


This subroutine sets up the matrices corresponding to tables I to IV. The assigned 
thermodynamic state being set up (tables I and II) is specified by the following codes: 


Assigned 

Codes 

thermodynamic 


state 


TP 

TP = .TRUE. VOL = .FALSE. CONVG = FALSE. 

HP 

HP = .TRUE. VOL = .FALSE. CONVG = .FALSE. 

SP 

SP = .TRUE. VOL = .FALSE. CONVG = .FALSE. 

TV 

TP = .TRUE. VOL = .TRUE. CONVG = FALSE. 

UV 

HP = .TRUE. VOL = .TRUE. CONVG = .FALSE. 

sv 

SP = .TRUE. VOL = .TRUE. CONVG = .FALSE. 


After convergence of any of the previous six problems, setup of the derivative ma 
trices (tables III and IV) is specified by the following codes: 


Derivative 

Codes 

DLVTP 

DLVPT 

CONVG = .TRUE. LOGV = .FALSE. 
CONVG = .TRUE. LOGV = .TRUE. 
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Subroutine GAUSS 


Subroutine GAUSS is used to solve the set of simultaneous linear iteration equations 
constructed by subroutine MATRIX. The solution is effected by performing a Gauss re- 
duction using a modified pivot technique. In this modified pivot technique only rows are 
interchanged. The row to be used for the elimination of a variable is selected on the 
basis that the largest of its elements, after division by the leading element, must be 
smaller than the largest element of the other rows after division by their leading ele- 
ments . 

The solution vector is stored in X(k). In the event of a singularity, IMAT (which is 
equal to the number of rows) is set equal to IMAT - 1. IMAT is tested later in subrou- 
tine EQLBRM. 


Subroutine 0UT1 

This subroutine, together with entries OUT2andOUT3, writes statements common 
to all problems. OUT1 writes statements giving the data on REACTANTS and on o/f, 
percent fuel, equivalence ratio, and density. 

Entry OUT2 . - This entry writes the statements for printing values of pressure, 
temperature, density, enthalpy, entropy, molecular weight, ( d In V/3 In P)rp (if equilib- 
rium), (3 In V/0 In T)p (if equilibrium), heat capacity, yg, and sonic velocity. These 
variables and corresponding labels are printed with a variable format described in 
BLOCK DATA. 

Entry OUT3 . - Entry OUT3 writes statements giving the equilibrium mole fractions 
of reaction species. 


Subroutine VARFMT 

Subroutine VARFMT (variable format) adjusts the number of decimal places printed 
in F-format in the variable format, FMT, according to the size of the number. It is 
used for P c /P e » p » and A g /A t . Variable format is described in BLOCK DATA. 


Subroutine EFMT 

Subroutine EFMT (E -format) writes statements in a special exponent form. This 
form is similar to the standard FORTRAN E -format, but the letter E and some of the 
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spaces have been removed for compactness. It is used to write density and mole frac- 
tions with the TRACE option. 


Subroutine THERMP 

This subroutine is the application module for TP, HP, SP, TV, UV, and SV prob- 
lems. Common variables which must be set according to the assigned thermodynamic 
states are given in the section Apolication Modules (p. 76). For these problems, the 
variables TP, HP, SP, SO, and VOL are set or read in in the main program. HSUBO is 
set either in NEWOF or HCALC. The general flow of the routine is given in figure 3. 

Indices run from 1 to NP (NP ^ 26) both for assigned pressures P and assigned 
volumes (V in INPT2 and vL in THERMP). Indices rim from 1 to NT (NT < 26) for as- 
signed temperature T. NP and NT are set in the main program. 


Subroutine ROCKET 

This subroutine is the control program for the RKT problem (rocket performance 
calculations discussed in section ROCKET PERFORMANCE). A flow diagram for this 
subroutine is given in figure 5. Subroutine ROCKET obtains the required thermodynamic 
properties for equilibrium performance by calling subroutine EQLBRM. For frozen per- 
formance, subroutine ROCKET calls subroutine FROZEN to obtain the required thermo- 
dynamic properties. Rocket performance parameters are then obtained by calling sub- 
routine RKTOUT . In addition to calling RKTOUT and FROZEN, and in addition to using 
controls common to all problems (discussed in section MODULAR FORM OF THE PRO- 
GRAM, p. 75), subroutine ROCKET also does the following: 

(1) It reads and processes the input data in RKTINP namelist. 

(2) It calculates estimates for throat pressure ratios. 

(3) It calculates estimates for pressure ratios corresponding to assigned area ratios 
(if any). 


aotcTNALi PAGE IS 

Subroutine RKTOUT “ poor QUaUT* 

This subroutine calculates various rocket performance parameters from previously 
calculated thermodynamic properties. 

It is also the control program for writing rocket performance output. It contains the 
WRITE statements that apply specifically to rocket parameters and it calls subroutine 
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OUT1 and entries OUT2 and OUT3 for the WRITE statements common to all problems. 
The rocket parameters are printed with the variable format, FMT, described in BLOCK 

DATA. 

Subroutine RKTOUT is called from subroutine ROCKET . 


Subroutine FROZEN 

Subroutine FROZEN is called from ROCKET to calculate the temperature and ther- 
modynamic properties for the following assigned conditions: 

(1) Composition frozen at combustion conditions (NFZ = 1) 

(2) An assigned exit pressure (PP) 

(3) An assigned entropy equal to the entropy at combustion conditions (SO = SSUM(l)) 
The iteration procedure used for obtaining the exit temperature is discussed in the sec - 
tion Procedure for Obtaining Frozen Rocket Performance (p. 40). 

If a temperature is reached 50 K below the range of a condensed combustion species 
(TEMP(j, 1) to TEMP(j, 2)), calculations are stopped. TT is set to zero and control is 
returned to ROCKET where a message is printed and data for all preceding points are 
listed. 

The variables which must be set in common before calling FROZEN include NFZ, 
NPT, TT, PP, IUSE(j) , COEF(k,i,j), SO, NS, NC, TEMP(j,k), TMID, WM(1), EN(j,l), 
RR, and TOTN(l). The variables which are set by FROZEN include TTT(NPT), 
PPP(NPT) , SSUM(NPT) , HSUM(NPT) , CPR(NPT), GAMMAS(NPT) , VLM(NPT), 
WM(NPT) ,' DLVPT(NPT) , DLVTP(NPT), and TOTN(NPT) . Many of these variables are 
the same as are required by or set by EQLBRM . 


Subroutine SHCK 

Subroutine SHCK is the application module for the SHOCK problems. It calculates 
the shock parameters discussed in the section "INCIDENT AND REFLECTED SHOCKS' 
It reads and processes the input data in SHKINP namelist. Depending on which options 
are specified, it calculates incident shock conditions based on compositions frozen at 
.initial conditions and/or based on equilibrium compositions after shock. It also calcu- 
lates, based on specified options, frozen and/or equilibrium reflected shock conditions 
relative to equilibrium and/or frozen incident shock conditions. 
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Subroutine DETON 


This subroutine does the calculations required to obtain Chapman- Jouguet detonation 
properties as described in the section CHAPMAN -JOUGUET DETONATIONS. Detonation 
calculations are limited to gaseous reactants only. When no T schedule is included in 
the namelist input INPT2, enthalpies for the reactants may be either punched on the 
cards or calculated by subroutine HCALC where the zeros are punched in card columns 
37 and 38 as usual. However, if a T schedule is used, all enthalpies will be calculated 
by the program regardless of the punches on the cards. In this case, of course, all re- 
actants must be in the THERMO data. 

BLOCK DATA 

BLOCK DATA contains atomic data stored in ATOM(i, j) and many of the variables 
used with the variable format, FMT. The ATOM variables are defined in appendix B. 
The format variables are stored in the common labeled OUPT and are described here. 

A variable format was used so that one format, FMT, could be used in the final out- 
put with changes in the number of decimal places according to the sizes of the numbers . 
The format is used to print a label and from 1 to 13 associated numbers. The labels con- 
tain 14 alphameric characters stored in four words and printed with 3A4,A2. The num- 
bers are all printed in a field of 9. FMT is initially set in BLOCK DATA as follows: 

FMT (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) 

(1H , 3A4 ,A2, F9. 0, F9. 0, F9. 0, F9. 0, F9. 0, F9. 0, 

FMT (16) (17) (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30) 

F9. 0, F9. 0, F9. 0, F9. 0, F9. 0, F9. 0, F9. 0 ) 

where the spaces are stored as blanks. 

Some variables set in BLOCK DATA to modify FMT are as follows: 

Variable: F0 FI F2 F3 F4 F5 FB FMT13 FMT9X FMTI9 
Storage: 0, 1, 2, 3, 4, 5, 13 9X, 19, 

The following is a list of variables used as labels and printed with 3A4,A2 in FMT: 
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Variable 

Stored label 

FP 

P, ATM 

FT 

T. DEG K 

FH 

H, CAL G 

FS 

S, CAL (G)(K) 

FM 

M, MOL WT 

FV 

(DLV/DLP)T 

FD 

(DLV /DLT)P 

FC 

CP, CAL/(G)(K) 

FG 

GAMMA (S) 

FL 

SON VEL, M/SEC 

FR1 

PC/P 

FC1 

CF 

FN 

MACH NUMBER 

FR 

CSTAR, FT /SEC 

FI 

ISP, LB-SEC/LB 

FA 

IVAC, LB-SEC LB 

FA1 , FA2 

AE AT 


MODIFICATION OF THE PROGRAM 
Changing Number of Possible Reaction Species 

All variables dimensioned for the number of possible reaction species are in the 
common labeled SPECES which takes two cards. These variables are dimensioned for 
150 species in the main program and subroutines SEARCH, SAVE, HCALC, MATRIX, 
CPHS, EQLBRM , OUT1, THERMP, ROCKET, RKTOUT, FROZEN, SHCK, and DETON. 
The two cards are the same in all of these routines. The 150 number may be either in- 
creased or decreased, to any number as required. No other program changes are nec- 
essary. 

Actually the number 150 is considerably larger than necessary for most chemical 
systems. Since each species requires 51 storages, an obvious way to reduce storage is 
to decrease the size of this number. This procedure might be desirable if a new appli- 
cation were added and additional storage were required. 


Eliminating an Application 


Any application module may be removed simply by removing the statement calling 
the controlling subroutine (THERMP, ROCKET, SHCK, or DETON) and removing the 



subroutine(s) in the application module. The calling statements are near the end of the 
main program (see MODULAR FORM OF THE PROGRAM, p. 75). 


Adding an Application 

An application may be added by doing the following: 

(1) Giving the type of problem a name and including the name in the main program in 
the following places: 

(a) Namelist INPT2 

(b) Logical statement 

(c) Statement setting it to false before the INPT2 read 

(2) Programming an application module as described in the section Application 
Module (p. 76). 

(3) Calling the module when the name variable is true after the INPT2 namelist data 
have been processed in the main program. 


Modifying Input and Output Modules 

The input modules may be modified as long as the common variables used by 
EQLBRM are properly set. These variables are indicated in the description of the 
modules . 

The listed output may be modified. Labeled common variables which are calculated 
by the equilibrium module and which are used in subroutine OUT2 and OUT3 include 
HSUM, SSUM, CPR, DLVTP, DLVPT, GAMMAS, PPP, TTT, VLM, TOTN, and EN. 
(Variables are defined in appendix B.) 

The units in the output may be changed in OUT 2 and RKTOUT. Either a variable 
format (FMT) or a special E-format (subroutine EFMT) is used for all output. The vari- 
able format (FMT) is described in the section "BLOCK DATA" and the adjustments in 
number of decimal places may be easily reset in the output routines. The parameter 
labels may be adjusted in the DATA statements in BLOCK DATA, the output routines, or 
in the application module itself. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 23, 1971, 

129-01. 
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APPENDIX A 


SYMBOLS 


area, 

ratio of nozzle exit area to throat area, eq. (106) 
velocity of sound, m/sec, eqs. (67) and (74) 
velocity of sound, m/sec, defined by eq. (85) 

stoichiometric coefficients, kg-atoms of element i per kg-mole of 
species j, (kg -atom) ^(kg -mole) . 

Stoichiometric coefficients, kg-atoms of element i per kg-mole of 
reactant j (oxidant if k = 1, fuel if k = 2), (kg -atom) i y/(kg- mole ) ^ 

least squares coefficients, eqs. (90), (91), and (92) 

kg-atoms of element i per kg of mixture, (kg-atom^Ag, eq. (7b) 

assigned kg-atoms of element i per kilogram of total oxidant (k = 1) or 
total fuel (k = 2), (kg-atom)^/kg^ k \ eq. (187) 

assigned kg-atoms of element i per kg of total reactant, (kg-atom^Ag, 
eq. (191) 

coefficient of thrust, eq. (105) 

standard state constant pressure specific heat for species or reactant j , 
J/(kg-mole).(K) 

standard state constant volume specific heat for species j, J /(kg- mole) ^(K) 


constant pressure specific heat of total reactant, kg-moleAg, eq. (206) 


( c p) s / r ' di 


dimensionless constant pressure specific heat for species j 


dimensionless constant volume specific heat for species j 


characteristic velocity, m/sec, eq. (103) 

constant pressure specific heat of mixture, J/(kg)(K), eqs. (49), (59), 
or (121) 
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constant volume specific heat of mixture, J/(kg)(K), eq. (70) or (124) 

Helmholtz free energy of mixture with constraints, defined by eq. (33), 
J/kg; or force N, eq. (98) 

percent of total fuel in total reactant by weight (or mass) 

Helmholtz free energy of mixture, J/kg, eqs. (30) and (31) 

fuel to air weight (or mass) ratio or fuel to oxidant weight (or mass) ratio 

Gibbs free energy of mixture with constraints, J/kg, defined by eq. (8) 

Pj/RT, dimensionless Gibbs free energy for species j 

Gibbs free energy of mixture, J/kg, eq. (5) 

conversion factor 

standard state enthalpy for species j, J/(kg-moie)j 

enthalpy of reactant j (k = 1 for oxidant, k = 2 for fuel), J/(kg-mole)( k \ 
eq. (192) ] 

heat of formation at temperature T, J/(kg-mole) 

h/RT, enthalpy of mixture, (kg-mole)/kg 

(h°)/r T , dimensionless enthalpy for species j 

enthalpy of total reactant, kg-moleAg, eq. (193) 

enthalpy of total oxidant (k = 1) or total fuel (k = 2), (kg-moleAg) k , 
eq. (192) 

enthalpy of mixture, JAg, eq. (14) 

enthalpy of total reactants, JAg, eq. (193) 

term defined by right side of eq. (156) 

term defined by right side of eq. (171) 

term defined by right side of eq. (131) 

specific impulse, N/(kg/sec) or m/sec, eq. (99) 

specific impulse with exit and ambient pressures equal, N/(kg/sec) or 
m/sec, eq. (100) 

vacuum specific impulse, N/(kg/sec) or m/sec, eq. (101) 
molecular weight of mixture, kg/kg-mole, eq. (2) 
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P* 
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r 
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atomic weight of chemical element i, (kg /kg -atom) ^ eq. (190) 
molecular weight of total reactant, kg kg-mole, eq. (197) 

molecular weight of total oxidant (k = 1) or total fuel (k = 2) , 
(kgAg-mole/ k \ eq. (196) 

molecular weight of reactant j (oxidant if k = 1, fuel if k = 2), 
(kgAg-mole)j k \ eq. (190) 

Mach number, eq. (102) 

moles of reactant j per kilogram of total reactant, (kg -mole) -Ag, 
eq. (205) 

mass flow rate, kg /sec, eq. (96) 

(k) 

kg -moles of reactant j (oxidant if k = 1, fuel if k = 2), (kg-mole). 

moles of mixture, kg-mole Ag, eq* ( 3 ) 

kg-moles of species j per kg of mixture, (kg-mole)jAg 

kg- moles of reactant j per kg of total oxidant (k = 1) or total fuel (k = 2) , 
(kg-mole) j k) /kg (k) , eq. (188) or (189) 

oxidant to fuel weight (or mass) ratio 
2 

pressure, N/m 

pressure, atm (see eqs. (11) and (17)) 

ratio of combustion pressure to exit pressure 

2 

assigned or initial pressure, N/m 

term defined by right side of eq. (155) 

term defined by right side of eq. (170) 

term defined by right side of eq. (130) 

universal gas constant, 8314.3 J/(kg-mole)(K) 

equivalence ratio, eq. (204) 

entropy of species j, J/(kg-mole).(K), eq. (17) 

standard state entropy for species j, J/(kg-mole).(K) 
s/R, entropy of mixture, (kg-mole) Ag 
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Sj/R, dimensionless entropy for species j 
entropy of total reactants, (kg -mole) /kg, eq. (207) 
entropy of mixture, J/(kg)(K), eq. (16) 
entropy of total reactants, J/(kg)(K), eq. (207) 
temperature, K 

standard state internal energy for species j, J/(kg-mole). 

internal energy of reactant j (oxidant if k = 1, fuel if k = 2), 
(J/kg-mole)p, eq. (194) 

(u°)/rT, dimensionless internal energy for species j 
internal energy of total reactant, (kg-mole)/kg, eq. (195) 

internal energy of total oxidant (k = 1) or total fuel (k = 2) , 

(kg -mole/kg ) ^ k \ eq. (194) 

velocity (in shock problems, velocity relative to incident or reflected 
shock front) , m/sec 

velocity at station 1 (in shock problems, velocity of unshocked gas rela- 
tive to incident shock front) , m/sec 

velocity at station 2 (in shock problems, velocity of incident -shocked gas 
relative to incident shock front), m/sec 

velocity of incident-shocked gas relative to reflected shock front, m/sec, 
eq. (128) 

velocity of reflected-shocked gas relative to reflected shock front, m/sec, 
eq. (129) 

internal energy of mixture, J/kg, eq. (38) 
internal energy of total reactant, J/kg, eq. (195) 
specific volume, m /kg 

positive oxidation state of total reactant, eq. (202) 
negative oxidation state of total reactant, eq. (203) 

positive oxidation state of total oxidant (k = 1) or total fuel (k = 2) , 
eq. (200) 





v: 

i 

V 


V 1 

v 2 


v 5 


W j k) 


W R 

W S 

y 

y s 

y S,T 

X 



r 

Po 


negative oxidation state of total oxidant (k = 1) or total fuel (k - 2) , 
eq. (201) 

positive oxidation state of chemical element i 
negative oxidation state of chemical element i 

actual velocity of shocked or unshocked gases in fixed coordinates, m/sec, 
eq. (125) 

actual velocity of unshocked gases in fixed coordinates (vj = 0), eq. (126) 
actual velocity of incident- shocked gases in fixed coordinates, m/sec, 
eq. (127) 

actual velocity of reflected-shocked gases in fixed coordinates (v & as- 
sumed to be zero), m/sec, eq. (129) 

(k) 

weight of reactant j (oxidant if k = 1, fuel if k = 2), kg^ , eq. (188) 

shock front velocity, m/sec, eq. (125) 
reflected shock front velocity, m/sec, eq. (128) 
incident shock front velocity, m/sec, eq. (126) 
ratio of specific heats, eq. (72) 
isentropic exponent, eqs. (71) and (73) 
special isentropic exponent defined by eq. (84) 
control factor defined by eq. (78) or (146) 

Lagrangian multiplier for chemical element i, J/(kg-atom) i 

control factor defined by eq. (76) 

control factor defined by eq. (77) 

control factor defined by eq. (144) 

control factor defined by eq. (145) 

chemical potential of species j, J/(kg-mole)j, eq. (6) 

standard state chemical potential for species j, J/(kg-mole). 

-X./RT, Lagrangian multiplier for chemical element i, kg-mole/ 

(kg -atom) i 

O 

density of mixture, kg/m , eq. (1) 

o 

density of total reactant, kg/m , eq. (199) 
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Subscripts: 

a 

c 

e 


_ « (k) 

density of total oxidant (k = 1) or total fuel (k = 2), (kg/m J ) , eq. (198) 

, o (k) 

density of reactant j (oxidant if k = 1, fuel if k - 2), (kg/m' J ). , 

eq. (198) J 

ambient 

combustion; condensed 
exit 


g gas 

k iteration k 


m melting temperature 

o total reactant; zero** 1 iteration 

t throat 

1,2,5 stations 

Superscripts: 

0 an assigned or initial condition 

k 1 if oxidant, 2 if fuel 

Indices: 

1 number of chemical elements (if ions are considered, number of chemical 

elements plus one) 

m number of possible gaseous species 

n number of possible species, gases and condensed 

NREAC number of reactants 
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APPENDIX B 


COMMON VARIABLES 


All common variables are given in this appendix except those with common label 
OUPT. These latter variables are discussed in section BLOCK DATA (p. 90). 


Variable 

Dimension 

Common 

Routines 

Description and comments 



label 

where used 


A 

15,150 

SPECES 

SEARCH 

MATRIX 

EQLBRM 

HCALC 

a.., stoichiometric coefficient of i 

1 3 fii 

element in j species. 

AC 

2 

MISC 

SAVE 

Saves JSOL (AC(1)) and JLIQ 
(AC(2)) . 

AEAT 

13 

PERF 

ROCKET 

In RKTOUT, area ratios A e /A t to 




RKTOUT 

be printed out. In SHCK, M2M1, 




SHCK 

molecular weight ratio. In DETON, 




DETON 

GM1, y for unburned gas. 

AM 

2 

MISC 

REACT 

AM(k) s M^ k \ molecular weight of 




HCALC 

total oxidant (k = 1) or total fuel 
(k = 2) . 

ANUM 

15,5 

MISC 

REACT 

Stoichiometric coefficients of reac- 




HCALC 

tants. First subscript refers to the 




OUT1 

reactant, the second, the element on 




SHCK 

the card. (See section Subroutine 
REACT). 

APP 

13 

PERF 

ROCKET 

P c /P e , ratio of chamber pressure 




RKTOUT 

to exit pressure. In SHCK, RRHO, 




SHCK 

a density ratio. 

AREA 

1 

PERF 

ROCKET 

Logical variable indicating (if true) 
area ratios have been assigned. 
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Variable 

Dimension 

Common 

Routines 

Description and comments 



label 

where used 


ATOM 

3,101 

MISC 

REACT 

For j th atom: 




BLOCK DATA 

ATOM (l,j) = atomic symbol 
(H, HE, . . .). 

ATOM (2,j) = atomic weight. 
ATOM (3,j) = atomic valence. 

AWT 

1 

PERF 

ROCKET 

nT/Pu at throat. AWT is used in 
obtaining area ratios, eq. (106). 

E0 

15 

MISC 

NEWOF 

b?, assigned kg -atoms of element i 




MATRIX 

per kg of total reactant. 

BOP 

15,2 

MISC 

Main 

B0P(i,l) h bp\ kg -atoms of i** 1 




REACT 

NEWOF 

element per kg of total oxidant. 
B0P(i,2) s bp', kg -atoms of i** 1 





element per kg of total fuel, 
eq. (187). 

CALCH 

1 

INDX 

Main 

Logical variable indicating (if true) 




SAVE 

zeros were punched in card columns 




HCALC 

37 and 38 of at least one 




THERMP 

ROCKET 

SHCK 

REACTANTS card. 

COEF 

2,7,150 

SPECES 

SEARCH 

Constants in empirical equations for 




CPHS 

thermodynamic functions. 
COEF(i,j,k): i = 1 for upper tem- 
perature interval and i = 2 for the 
lower one; j = 1, . . . 7, for the 
7 coefficients, a.; k = 1, . . . 

150 for the number of species, eqs. 
(90) to (92) . 

CONVG 

1 

INDX 

EQLBRM 

Logical variable indicating conver- 




MATRIX 

gence is true. 


FROZEN 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

CPCVEQ 

1 

INDX 

EQLBRM 

MATRIX 

Equivalenced to LOGV in EQLBRM 
and MATRIX (see LOGV) 

CPCVFR 

1 

INDX 


Temporary storage. 

CPR 

13 

POINTS 

EQLBRM 

OUT2 

DETON 

SHCK 

ROCKET 

FROZEN 

c /R, constant pressure specific 
heat, kg-moleAg (eq. (59) or (121)). 

CPSUM 

1 

MISC 

EQLBRM 

CPHS 

MATRIX 

HCALC 

ROCKET 

FROZEN 

SHCK 

c /R, frozen constant pressure spe- 
cific heat, kg-moleAg, eq. (121). 

CR 

1 

MISC 


Not used. 

CSTR 

1 

PERF 

RKTOUT 

c* , characteristic velocity, ft/sec. 

DATA 

22 

MISC 

Main 

HCALC 

REACT 

SHCK 

DETON 

Temporary storage. 

DELN 

150 

SPECES 

EQLBRM 

Main 

In EQLBRM: Corrections to In n. 
for gases or n^ for condensed. 

In Main: Equivalent to ENSERT, an 
array of names read from INSERT 
cards. 

DENS 

15 

MISC 

REACT 

OUT1 

Density of reactant read from input 
card. 
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Variable 

DLVPT 

DLVTP 

EN 


ENLN 

ENLSAV 

ENN 


Dimension Common Routines Description and comments 

label where used 


13 POINTS EQLBRM 

OUT2 

DETON 

SHCK 

FROZEN 

13 POINTS EQLBRM 

OUT2 

DETON 

SHCK 

ROCKET 

FROZEN 

150,13 SPECES Main 

EQLBRM 

CPHS 

MATRIX 

SAVE 

OUT3 

SEARCH 

SHCK 

HCALC 

RKTOUT 

FROZEN 



eq. (51) or (123). 


eq. (50) or (122). 


i 


I 


nj, kg-moles of j* speciesAg of 
mixture (j < 150). The second sub- 
script is for point number for out- 
put. In SEARCH, EN is equiva- 
lenced to DATE, the dates read 
from THERMO data. 


150 


SPECES Main 

SEARCH 

EQLBRM 

MATRIX 

SAVE 


In n. (see EN). In the main pro- 
gram ENLN is equivalenced to 
OMIT, an array for storing names 
on OMIT cards. 


1 MISC SAVE 

1 MISC Main 

EQLBRM 

MATRIX 

SAVE 

ROCKET 


In n, saved for a particular point 
(see ENN). 

n, moles of mixture, kg-moleAg 
(eq. (3) for assigned volume prob- 
lems or eq. (79) for assigned pres 
sure problems) . 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

ENNL 

1 

MISC 

Main 

EQLBRM 

SAVE 

In n, (see ENN). 

ENSAVE 

1 

MISC 

SAVE 

n, saved for a particular point (see 
ENN). 

ENTH 

15 

MISC 

REACT 

HCALC 

OUT1 

For assigned pressure problems, 
(h£) , assigned enthalpy of a reac- 
tant in calories per mole. For as- 
signed volume (or density) prob- 
lems, , assigned internal 

energy of a reactant in calories per 
mole. 

EQL 

1 

PERF 

Main 

ROCKET 

OUT3 

RKTOUT 

SHCK 

DETON 

Logical variable indicating (if true) 
equilibrium calculations rather than 
frozen. 

EQRAT 

1 

MISC 

Main 

NEWOF 

OUT1 

r, equivalence ratio (eq. (204)). 

FAZ 

15 

MISC 

Main 

REACT 

HCALC 

OUT1 

Alphameric letter read from reac- 
tant cards. The letter G indicates 
gas; any other letter indicates con- 
densed phase . 

FOX 

15 

MISC 

REACT 

HCALC 

OUT1 

Alphameric letter, F or O, read 
from reactant cards indicating reac- 
tant is a fuel (F or blank) or oxi- 
dant (O) . 

FPCT 

1 

MISC 

Main 

Logical variable in namelist INPT2 
(see table VTII) . 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

FROZ 

1 

PERF 

ROCKET 

Logical variable in namelist 
RKTINP. If true, rocket perform- 
ance is desired with frozen compo- 
sition. 

G 

20,21 

(Double 

precision) 

DOUBLE 

EQLBRM 

MATRIX 

GAUSS 

SHCK 

Matrix region. 

GAMMAS 

13 

POINTS 

EQLBRM 

OUT2 

DETON 

SHCK 

ROCKET 

RKTOUT 

FROZEN 

y s or yg j isentropic exponent ex- 
ponent (eq. (71), (73), or (84)). 

HO 

150 

SPECES 

EQLBRM 

CPHS 

MATRIX 

HCALC 

FROZEN 

RKTOUT 

SHCK 

OUT3 

Jfj, dimensionless enthalpy for 
species j. HO is used as a tempor- 
ary storage in OUT 3. 

HP 

1 

INDX 

Main 

EQLBRM 

MATRIX 

THERMP 

DETON 

ROCKET 

Logical variable indicating (if true) 
either enthalpy and pressure or in- 
ternal energy and volume (or den- 
sity) have been assigned. 

HPP 

2 

MISC 

REACT 

HCALC 

NEWOF 

or ♦^T, enthalpy or inter- 
nal energy of total oxidant (k = 1) and 
total fuel (k = 2) (eqs. (192) and 
(194)). 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

HSUBO 

1 

MISC 

NEWOF 

HCALC 

MATRIX 

DETON 

SHCK 

THERMP 

h^R or u^/R assigned enthalpy or 
zero internal energy of total reac- 
tant (eqs. (193) and (195)). 

HSUM 

13 

POINTS 

EQLBRM 

MATRIX 

OUT2 

DETON 

SHCK 

ROCKET 

RKTOUT 

FROZEN 

in MATRIX and h/R at the end 
of EQLBRM and in the other rout- 
tines (eqs. (29) and (14)) . 

IC 

1 

INDX 


Not used. 

IDEBUG 

1 

INDX 

Main 

EQLBRM 

DETON 

SHCK 

ROCKET 

THERMP 

When set to a value i (i =» 0), inter- 
mediate output is printed for all 
points 5: i. 

IMAT 

1 

INDX 

EQLBRM 

MATRIX 

GAUSS 

Number of rows in matrix setup. 

IONS 

1 

INDX 

Main 

EQLBRM 

SAVE 

Logical variable indicating (if true) 
ionic species are to be considered. 

IP 

1 

INDX 

THERMP 

DETON 

ROCKET 

Index for pressure or specific 
volume values. 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

IQ1 

1 

INDX 

Main 

EQLBRM 

SAVE 

MATRIX 

Number of elements plus the number 
of condensed species included in ma- 
trix equations plus 1 . 

IQSAVE 

1 

INDX 

SAVE 

IQ1 saved for a particular point 
when ISV is negative . 

ISUB 

1 

INDX 

ROCKET 

Index for SUBAR values (see 
SUBAR). 

ISUP 

1 

INDX 

ROCKET 

Index for SUPAR values (see 
SUPAR) . 

ISV 

1 

INDX 

SAVE 

THERMP 

DETON 

SHCK 

ROCKET 

ISV > 0: Use compositions from 
last point (numbered ISV) as esti- 
mates for current point (numbered 
NPT). ISV < 0: Store compositions 
in SLN array to use as estimates for 
future points. Make ISV positive and 
proceed as with ISV > 0. ISV = 0: 
Use compositions stored in SLN 
array as estimates for current point 
(numbered NPT) . 

IT 

1 

INDX 

Main 

THERMP 

ROCKET 

RKTOUT 

DETON 

Index for T values (see T) . 

ITM 

1 

INDX 

ROCKET 

FROZEN 

Integer 1. 

ITNUM 

1 

INDX 

ROCKET 

Iteration number for converging to 
assigned area ratio. 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

IUSE 

150 

SPECES 

Main 

EQLBRM 

CPHS 

MATRIX 

HCALC 

SAVE 

SEARCH 

FROZEN 

SHCK 

See section "IUSE array”. 

JLIQ 

1 

INDX 

EQLBRM 

SAVE 

Index of condensed species which is 
included simultaneously with another 
condensed phase of the same spe- 
cies. JLIQ is the index for the 
phase in the higher temperature 
range; JSOL is the index for the 
phase in the lower temperature 
range. Otherwise JLIQ = JSOL = 0. 

JS1 

1 

INDX 

CPHS 

EQLBRM 

FROZEN 

HCALC 

SHCK 

Index for species for which ther- 
modynamic functions are to be 
calculated in CPHS . 

JSOL 

1 

INDX 

EQLBRM 

ROCKET 

SAVE 

See JLIQ. 

KASE 

1 

INDX 

Main 

OUT1 

Optional identifying case number for 
a set of reactants. The variable is 
part of namelist INPT2. 

KMAT 

1 

INDX 

EQLBRM 

MATRIX 

Number of columns in matrix setup. 
KMAT = IMAT + 1 . 
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Variable 

Dimension 

Common 

Routines 

Description and comments 



label 

where used 


LLMT 

15 

MISC 

Main 

Alphameric symbols for elements 




REACT 

in the chemical system being 




SEARCH 

HCALC 

EQLBRM 

SAVE 

NEWOF 

processed. 

LOGV 

1 

INDX 

EQLBRM 

Logical variable indicating (if true) 




MATRIX 

matrix setup given in table IV . 
LOGV is equivalenced to CPCVEQ 
which appears in common INDX. 

LSAVE 

1 

INDX 

SAVE 

LLMT(NLM), the last element sym- 
bol is saved in LSAVE for use at 
some future point. If IONS is true 
and LSAVE is not "E", then ionic 
species were removed. 

MOLES 

1 

INDX 

Main 

Logical variable indicating (if true) 




REACT 

that the relative amounts of the 




OUT1 

reactants have been specified in 




HCALC 

terms of moles. 

NAME 

15,5 

MISC 

Main 

Alphameric symbols for elements as 




REACT 

given on the reactant cards. Allow- 




HCALC 

ance is made for up to 5 symbols 




DETON 

SHCK 

OUT1 

and 15 reactants. 

NC 

1 

INDX 

Main 

Number of condensed species in 




SEARCH 

thermodynamic data for a particu- 




EQLBRM 

FROZEN 

lar system. 
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Variable 

NEWR 

NFZ 

NLM 

NOF 

NOMIT 

NP 

NPP 


Dimension Common 
label 


Routines 
where used 


1 


INDX Main 

REACT 

SEARCH 


1 


1 


1 


1 


1 


1 


INDX ROCKET 

FROZEN 
RKTOUT 

INDX Main 

REACT 

SEARCH 

HCALC 

NEWOF 

EQLBRM 

MATRIX 

INDX Main 

THERMP 

DETON 

SHCK 

ROCKET 

INDX Main 

SEARCH 

EQLBRM 

MATRIX 

REACT 

INDX Main 

THERMP 

DETON 

ROCKET 

INDX ROCKET 


Description and comments 


Logical variable indicating (if true) 
that the search of tape 4 has not yet 
been made for thermodynamic data 
for species in a new chemical sys- 
tem. 

Index of point where equilibrium 
composition is to be frozen. 

NFZ = 1 for rocket problems. 

Number of elements in the chemical 
system. 


Number of o/f values. Note: 
mixture values given in the form of 
equivalence ratio, percent fuel, or 
fuel-air ratio are converted to o/f 
values in the main program. 

Number of species to be omitted. 


Number of pressure values in the P 
array in namelist INPT2. 


Number of P /P„ values. This is 
c e 

equal to the number of PCP values 
in namelist RKTINP plus 2 to in- 
clude chamber and throat. 
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Variable 

NPT 

NREAC 

NS 


NSERT 
NSUB 
NS UP 
NT 


Dimension 

Common 

label 

Routines 
where used 

Description and comments 

1 

INDX 

All routines 
except Main, 
REACT, 
GAUSS, and 
SEARCH 

Index for the data saved for output 
(1 <; NPT < 13). 

1 

INDX 

Main 

REACT 

HCALC 

OUT1 

DETON 

SHCK 

Number of reactants. 

1 

INDX 

Main 

SEARCH 

HCALC 

OUT3 

EQLBRM 

CPHS 

MATRIX 

FROZEN 

SHCK 

RKTOUT 

SAVE 

Number of species considered in a 
particular chemical system. 

1 

INDX 

Main 

Number of condensed species listed 
on INSERT cards. 

1 

INDX 

ROCKET 

Number of SUBAR values (see 
SUBAR). 

1 

INDX 

ROCKET 

Number of SUPAR values (see 
SUPAR). 

1 

INDX 

Main 

THERMP 

DETON 

ROCKET 

Number of temperature values 
listed in T schedule in namelist 
INPT2. 
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Variable 

OF 

OXF 

P 

PCP 

PECWT 

PP 


Dimension 

Common 

label 

Routines 
where used 

Description and comments 

1 

MISC 

Main 

HCALC 

NEWOF 

OUT1 

DETON 

SHCK 

ROCKET 

THERMP 

o/f, oxidant to fuel weight ratio for 
current point. 

15 

MISC 

Main 

THERMP 

DETON 

SHCK 

ROCKET 

o/f, oxidant to fuel weight ratios. 

In the main program, OXF is equiv- 
alenced to MIX, an ENPT2 namelist 
variable . 

26 

POINTS 

Main 

THERMP 

DETON 

SHCK 

ROCKET 

Assigned pressure schedule in 
namelist INPT2. Values are con- 
verted to atmosphere units in the 
main program. For TV, UV, and 
SV problems, assigned volumes, 
VL, in cmVg are stored in the P 
array in the main program. 

22 

PERF 

ROCKET 

SHCK 

DETON 

Pressure ratios read in with 
RKTINP namelist. Storage is used 
as temporary storage in SHCK and 
DETON. 

15 

MISC 

REACT 

HCALC 

OUT1 

Relative weight or number of moles 
of reactants as read from card col- 
umns 46 to 52 on the reactant cards 

1 

MISC 

HCALC 

EQLBRM 

THERMP 

DETON 

SHCK 

ROCKET 

FROZEN 

Assigned pressure in atmospheres 
for current point. 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

PPP 

13 

POINTS 

EQLBRM 

OUT2 

DETON 

SHCK 

RKTOUT 

FROZEN 

Static pressure in atmospheres 
stored for output. 

R 

1 

MISC 

Main 

REACT 

HCALC 

OUT2 

THERMP 

DETON 

Universal gas constant, 1.987165 
cal/(g-mole)(K). 

RH 

2 

MISC 

REACT 

NEWOF 

p^\ density of total oxidant and 
(2) 

p , density of total fuel, eq. (198) 

RHOP 

1 

MISC 

Main 

OUT1 

NEWOF 

p Q , density of total reactant, 
eq. (199). 

RMW 

15 

MISC 

REACT 

HCALC 

Molecular weight of individual 
reactants. 

RR 

1 

MISC 

Main 

OUT2 

EQLBRM 

DETON 

SHCK 

ROCKET 

RKTOUT 

FROZEN 

Universal gas constant, 8314.3 
J/(kg-mole)(K). 

RTEMP 

15 

MISC 

REACT 

HCALC 

OUT1 

DETON 

SHCK 

Temperatures read from reactant 
card columns 64 to 71. 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

S 

150 

SPECES 

HCALC 

EQLBRM 

CPHS 

MATRIX 

FROZEN 

SHCK 

(3°) Ir, dimensionless entropy 

calculated from coefficients for 
reaction species. See appendix A 
for definition. 

SHOCK 

1 

INDX 

Main 

HCALC 

EQLBRM 

Logical variable indicating (if true) 
SHOCK problem. 

SLN 

150 

SPECES 

SAVE 

In n. for gases and n^ for con- 
densed species for a particular 
point (see ISV with negative value) . 

SO 

1 

MISC 

Main 

EQLBRM 

MATRIX 

ROCKET 

FROZEN 

In INPT2 namelist for SP and SV 
problem, entropy in cal/(g)(K). For 
SHOCK problems SO = J Q , 
eq. (207). For RKT problems, 

SO = s/R, eq. (16), for the 
combustion point (NPT = 1) . 

SONVEL 

13 

POINTS 

OUT2 

DETON 

RKTOUT 

a, velocity of sound, eq. (74), 
or ag rp, eq. (85). 

SP 

1 

INDX 

Main 

EQLBRM 

MATRIX 

THERMP 

ROCKET 

Logical variable indicating (if true) 
entropy and either pressure (VOL 
is false) or volume (VOL is true) 
ha vp been assigned. In namelist 
INPT2, SP indicates an SP problem 

SPIM 

13 

PERF 

DETON 

SHCK 

RKTOUT 

In RKTOUT, I , specific impulse 
with ambient and exit pressures 
equal, (lb-force) (sec) /lb-mass. 
Storage is used for other variables 
in SHCK and DETON. 

SSO 

1 

PERF 


Not used. 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

SSUM 

13 

POINTS 

HCALC 

EQLBRM 

OUT2 

ROCKET 

FROZEN 

SUCK 

s/R, entropy of mixture, 

(kg -mole) Ag- 

SUB 

150,3 

SPECES 

Main 

SEARCH 

EQLBRM 

OUT1 

RKTGUT 

SHCK 

Alphameric names of species read 
from THERMO data. Each name 
takes 3 words, 4 characters each. 

SUBAR 

13 

PERF 

DETON 

SHCK 

ROCKET 

In ROCKET and RKTINP namelist, 
subsonic area ratios. SUBAR is 
equivalenced to other variables in 
SHCK and DETON. 

SUMN 

1 

MISC 

Main 

EQLBRM 

MATRIX 

m 

nj, moles of mixture, 

kg -mole Ag- 

SUPAR 

13 

PERF 

DETON 

SHCK 

ROCKET 

in ROCKET and RKTINP namelist, 
supersonic area ratios. SUPAR is 
equivalenced to other variables in 
SHCK and DETON. 

T 

26 

POINTS 

Main 

THERMP 

DETON 

SHCK 

ROCKET 

RKTOUT 

T, assigned temperature schedule 
in namelist INPT2. Values are in 
degrees Kelvin. 

TEMP 

50,2 

SPECES 

SEARCH 

HCALC 

EQLBRM 

GAUSS 

FROZEN 

Temperature ranges of individual 
condensed species. The lower value 
is in TEMP(j, 1) . 
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Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

THIGH 

1 

MISC 

Main 

SEARCH 

HCALC 

EQLBRM 

SHCK 

Highest temperature in temperature 
range of species in the THERMO 
data. 

TLN 

1 

MISC 

EQLBRM 

CPHS 

HCALC 

SHCK 

FROZEN 

In T. 

TL/OW 

1 

MISC 

Main 

SEARCH 

HCALC 

EQLBRM 

FROZEN 

SHCK 

Lowest temperature in temperature 
range of species in the THERMO 
data. 

TM 

1 

MISC 

HCALC 

EQLBRM 

MATRIX 

ln<P atm /n) ' 

TMID 

1 

MISC 

Main 

SEARCH 

CPHS 

Common temperature of two temper- 
ature ranges of the THERMO data. 

TOTN 

13 

POINTS 

EQLBRM 
OUT 3 
FROZEN 
RKTOUT 

NS 

iij, used in the calculation of 

i=i 

mole fractions. 

TP 

1 

INDX 

Main 

EQLBRM 

MATRIX 

THERMP 

DETON 

SHCK 

ROCKET 

Logical variable indicating (if true) 
temperature and either pressure 
(VOL is false) or volume (VOL is 
true) have been assigned . The as- 
signments may be for a complete 
problem (TP) or a particular point. 
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Variable 

TRACE 

TT 


TTT 


Dimension Common Routines 


label where used 


Description and comments 


1 


1 


MISC Main 

EQLBRM 
OUT 3 
RKTOUT 
SHCK 

MISC NEWOF 

HCALC 
EQLBRM 
CPHS 
MATRIX 
THERMP 
DETON 
SHCK 
ROCKET 
FROZEN 


If TRACE >0, print mole frac- 
tions 2 : TRACE in special E -format. 
TRACE is namelist INPT2 variable. 


Current temperature . 


13 


13 


POINTS EQLBRM 
OUT 2 
DETON 
SHCK 
ROCKET 
RKTOUT 
FROZEN 

POINTS Main 

REACT 

OUT2 

OUT3 

VARFMT 

THERMP 

DETON 

SHCK 

RKTOUT 

EFMT 


Temperature stored for printout. 


Temporary storage. Assigned 
volume in INPT2 namelist. 
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Variable 

VACI 

VLM 

VMIN 

VMOC 

VOL 

VPLS 

WM 


Dimension Common Routines Description and comments 

label where used 


13 


13 


2 


PERF RKTOUT 
DETON 
SHCK 


POINTS EQLBRM 
THERMP 
SHCK 
FROZEN 
OUT2 

MISC Main 

REACT 

NEWOF 


In RICTCUT , I vac , (lb force) 
(sec)/(lb mass) Stor age is used 
for other variables in SUCK and 
DETCN. 

Volume, cm 3 /g, stored for output. 


V"^ , negative oxidation state of 
total oxidant (if k = 1) or total fuel 
(if k = 2), eq. (201). 


13 

1 


2 


PERF 

RKTOUT 

uf , Mach number, eq. (102). Used 


DETON 

as temporary storage in DETON . 

INDX 

Main 

Logical variable indicating (if true) 


HCALC 

NEWOF 

EQLBRM 

MATRIX 

THERMP 

volume has been assigned. 

MISC 

Main 

V + ^ k \ positive oxidation state of 


REACT 

total oxidant (if k = 1) or total fuel 


NEWOF 

(if k = 2), eq. (200). 


13 


POINTS OUT2 

EQLBRM 

DETON 

SHCK 

ROCKET 

RKTOUT 

FROZEN 


M, molecular weight of the mixture, 
kgAg-mole, eq. (2). 


ORIGINAL PAGE IS 
OF POOR QUALITY 


117 



Variable 

Dimension 

Common 

label 

Routines 
where used 

Description and comments 

WP 

2 

MISC 

Main 

REACT 

HCALC 

NREAC 

WP(k) = y w( k) . 
j=l 

X 

20 (Double 
precision) 

DOUBLE 

EQLBRM 

GAUSS 

MATRIX 

SHCK 

Answer region for matrix solution 
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APPENDIX D 


THERMO DATA (FORMAT AND LISTING) 


The order and format of the input data cards in this appendix are given in the follow- 
ing table: 


Card 

order 

Contents 

Format 

Card 

column 

1 

THERMO 

3A4 

1 to 6 

2 

Temperature ranges for 2 sets of coefficients: 
lowest T, common T, and highest T 

3F10.3 

1 to 30 

3 

Species name 

3A4 

1 to 12 


Date 

2A3 

19 to 24 


Atomic symbols and formula 

4(A2, F3. 0) 

25 to 44 


Phase of species (S, L, or G for solid, liquid, 
or gas, respectively) 

A1 

45 


Temperature range 

2F10.3 

46 to 65 


Integer 1 

115 

80 

4 

Coefficients a^i = 1 to 5) in equations (90) to (92) 
(for upper temperature interval) 

5(E15. 8) 

1 to 75 


Integer 2 

15 

80 

5 

Coefficients in equations (90) to (92) (a g , for 

upper temperature interval, and a^ , a 2 , and 
for lower) 

5(E15. 8) 

1 to 75 


Integer 3 

15 

80 

6 

Coefficients in equations (90) to (92) (a., a-, a c 
3 .,. for lower temperature interval) 

4(E15. 8) 

1 to 60 

(a) 

Integer 4 

Repeat cards numbered 1 to 4 in cc 80 for each 
species 

120 

80 

(Final 

card) 

END (Indicates end of thermodynamic data) 

3A4 

1 to 3 


Gaseous species and condensed species with only one condensed phase can be in 
any order. However, the sets for two or more condensed phases of the same 
species must be adjacent. If there are more than two condensed phases of a 
species, their sets must be either in increasing or decreasing order according 
to their temperature intervals. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 



The program writes the THERMO data on tape 4 as the cards are read. The data are 
written on tape with the same formats as the cards with two minor exceptions. The 
THERMO code card is omitted and the card numbers in card colunn 80 are also 
omitted . 

A listing of the THERMO data cards included with the program follows . These cards 
were made by the program described in reference 15. 
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TABLE I. - ITERATION EQUATIONS TO DETERMINE EQUILIBRIUM COMPOSITIONS FOR EITHER ASSIGNED TEMPERATURE AND 






PRESSURE, ENTHALPY AND PRESSURE, OR ENTROPY AND PRESSURE 



^Column not uMd for assigned temperature and pressure. 
* > Row used only for assigned enthalpy and pressure. 
c Row used only for assigned entropy and pressure. 



TABLE II. - ITERATION EQUATIONS TO DETERMINE EQUILIBRIUM COMPOSITIONS FOR EITHER ASSIGNED TEMPERATURE AND 
VOLUME, INTERNAL ENERGY AND VOLUME, OR ENTROPY AND VOLUME 


(45) m m 
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Column not used for assigned temperature and volume. 
Row used only for assigned internal energy and volume 
Row used only for assigned entropy and volume. 
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TABLE III. - EQUATIONS FOR EVALUATING DERIVATIVES WITH RESPECT TO LOGARITHM 


OF TEMPERATURE AT CONSTANT PRESSURE 


^v^ariables 3 v j 

Equation'v. 3 In T 
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Equation 


TABLE IV. - EQUATIONS FOR EVALUATING DERIVATIVES WITH RESPECT TO 
LOGARITHM OF PRESSURE AT CONSTANT TEMPERATURE 


.^Variables I 3 
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m 
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- £ a 
j=i 


3 In n Right side 
a In P I 3 ln p 
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TABLE V - PROGRAM INPUT 


THERMO 
code card 


r °r 


THERMO 
data cards 



/reactants 

code card 


/ REACTANTS 
cards 


OMIT 

card{s) 


(optional) 


/ INSERT 
card(s) 


(optional) 


NAMELISTS 
code card 


$INPT2 

Optional variables: 

KASE 

MIX (1 to 15) 

OF, FPCT. FA, or ERATIO 

IONS 

IDEBUG 

TRACE 

(additional namelist input is 
given in table on the right) 


Problem 

Namelist 

Variables 



Required 

Optional 

Assigned temperature and 
pressure (TP) 

INPT2 

TP = .TRUE. 
T(1 to 26) 
P(1 to 26) 

NSQM, PS LA , or 
MMHG 

Assigned enthalpy and 
pressure (HP) 

INPT2 

HP = .TRUE. 
P(1 to 26) 

NSQM, PSIA, or 
MMHG 

Assigned entropy and 
pressure (SP) 

INPT2 

SP = .TRUE. 
S0(1) 

P(1 to 26) 

NSQM, PSIA, or 
MMHG 

Assigned temperature and 
volume or density (TV) 

INPT2 

TV = .TRUE. 
T(1 to 26) 
V(1 to 26) or 
RHO(l to 26) 


Assigned internal energy 
and volume or density (UV) 

INPT2 

UV = .TRUE. 
V(1 to 26) or 
RHO(l to 26) 


Assigned entropy and volume 
or density (SV) 

INPT2 

SV = .TRUE . 
V(1 to 26) or 
RHO(l to 26) 
S0<1) 


Detonation (DETN) 

INPT2 

DETN = .TRUE. 
P(1 to 26) finitiaf 

L * as J 

NSQM, PSIA. or 
MMHG 

T(1 to 26£nitial gas] 

Shock (SHOCK) 

INPT2 

SHOCK = .TRUE. 
P(1 to 13) [initial 
T(1 to 13) [ gas 

NSQM, PSIA, or 
MMHG 


SHKINP 

Ul(l to 13) or 
MACH1 (1 to 13) 

INCDEQ - .FALSE, 
or 

INCDFZ = . FALSE. 

Rocket (RKT) 

INPT2 

RKT = .TRUE. 
P(1 to 26) (cham- 
ber pressures) 

T(1 to 26)(cham- 
ber) NSQM, PSIA, 
or MMHG 


RKTINP 


EQL = . FALSE, or 
FROZ = . FALSE. 
PCP(1 to 22) 
SUPARO to 13) 
SUBARU to 13) 
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TABLE VI. - REACTANTS CARDS 


Order 

Contents 

Format 

Card columns 

First 

REACTANTS 

3A4 

1 to 9 

Any 

One card for each reactant species 
(maximum 15). Each card contains: 




(1) Atomic symbols and formula num- 
bers (maximum 5 sets) a 

5(A2, F7.5) 

1 to 45 


(2) Relative weight** or number of 
moles 

F7.5 

46 to 52 


(3) Blank if (2) is relative weight or 
M if (2) is number of moles 

A1 

53 


(4) Enthalpy or internal energy a , 
cal/mole 

F9.5 

54 to 62 


(5) State: S, L, or G for solid , 
liquid or gas, respectively 

A1 

63 


(6) Temperature associated with 
enthalpy in (4) 

F7.0 

64 to 70 


(6a) J If (4) is in units of kJ/kg-mole 
and blank if (4) is in units of 
cal/g-mole 

Al 

71 


(7) F if fuel or O if oxidant 

A1 

72 


3 

(8) Density in g/cm (optional) 

F8.5 

73 to 80 

Last 

Blank 




a Program will calculate the enthalpy or internal energy (4) for species in 
the THERMO data at the temperature (6) if zeros are punched in card 
columns 37 and 38. (See section Reactant enthalpy for additional in- 
formation.) 

^Relative weight of fuel in total fuels or oxidant in total oxidants. All 
reactants must be given either all in relative weights or all in number 
of moles. 




TABLE VO. - LIST OF REACTANTS CARDS FOR SOME OXIDANTS AND FUELS 
Chemical formula (card columns 1 to 45) Percent Assigned (a) Temper- (b) De 


kcc 46-92) enthalpy, 
cal /mole 


Acetonitrile 

Acetylene 

Air 

Aluminum 

Ammonia(g) 

Ammon ia(l) 

Ammonium perchlorate 

Aniline 

Argon 

Benzene 

Beryllium 

Butane 

1 -butene 

Chlorine (g) 

Chlorine (1) 

Chlorine trifluoride(g) 

Chlorine trifluorkle(l) 
Cyanogen(g) 

CyanogenO) 

Dlborane 
Bthnne 
Ethyl alcohol 
Ethylene 
Ethylene aside 
Ethylene polymer 
Fluor ine(g) 

Fluor lne(l) 

Graphite 

Helium 

Heptane 

Rydraalne 

*Ptaac: B, solid; L, liquid; C 
Sr net, F; ortdhnt, O. 
c Baaed on the following molar | 
*Estlmatc baaed an paraffin hy< 


N 1.S6176T .4153«J 4S.CC9324C .C0C1C0 
at l * 


ature, 

K 

(ce 64-71) 

t 29*. is 
L 192.60 
C 29®. 1 *» 

5 290. IS 

6 290.11 

t 2 J9. 1/ 

S 290. 1 1 
l 29®. IS 
G 29®. IS 
l 290. IS 
S 290. IS 
L 272.65 


(b) Density, 
g 'em 3 


References for References for 
assigned enthalpy density 


F .7057 
f .61 r 


r .0Tb 
f 1.95 
F 1.02173 


19 

14 (aeries E) 

(0 J 

reference element 
6 

14 (aeries E) 
14 (aeries A) 
21 

reference element 
22 

reference element 
22 


20 

14 (series E) 


14 (series E) 
14 (series A) 
20 


ICO. 

•5003. 

L 

266.9 

F 

.0243 

22 

22 

ICO. 

0. 

G 

290.15 

C 


reference element 


100* 

-5391. 

L 

239.39 

0 

1.54 

14 (serie® E) 

14 (scries E) 

ICO. 

-390C0. 

G 

290.15 

c 


19 


100. 

-45600. 

L 

2*4. 5*> 

0 

1.0517 

23 

24 

ICO. 

♦ 73040* 

G 

290.15 

F 


19 


100. 

67655. 

L 

252.01 

F 

.9517 

1, 15, 19, 25 

16 

100. 

497C » 

L 

10C.5 » 

F 

.4371 

14 (aerie® E) 

14 (aeriea E) 

100. 

-250C®. 

L 

104.5? 

F 

.5464 

21, 22 

22 

too. 

-66373. 

L 

296.15 

f 

.7093 

19 

SO 

100. 

0100. 

L 

169.44 

f 

.5600 

14 (aeries E) 

14 (series E) 

100. 

-10040. 

L 

203.72 

F 

• 0026 

19, 91, 95 

20 

100. 

-6100. 


290. IS 

F 

.935 

(d) 

27 

ICO. 

0. 

fr 

290.15 

0 


reference element 


too. 

-309®. 

L 

15. 32 

0 

1.505 

14 (aeriee D) 

14 (aeries D) 

too. 

0. 

5 

290.15 

F 

2.25 

reference element 

20 

ICO. 

0. 

0 

290.1s 

F 


reference element 


100. 

-43630. 

L 

290.15 

f 

.67951 

22 

22 

100* 

12100. 

t 

290.15 

f 

1.0036 

19 

21 


percents: N, - II. MSI, O, • 20. MM, A t ■ 1.1224, CO, • 0.0100. 
ttroearbon aeries. 




233 




TABLE VII - Concluded LIST OF REACTANTS CARDS FOR SOME OXIDANTS AND FUELS 


Chemical 

Chemical formula (card columns 1 to 45) 

Psrcsnt 
<cc 46-52) 

Assigned 
enthalpy, 
cal /mole 
(cc 54-62) 

(a) 

Temper - 
attire, 

K 

(cc 6-1-71) 


Density, 
g 'cm 3 
(cc 73-00) 

References for 
assigned enthalpy 

References for 
density 

Hydrogen<g) 

H 2. 


too. 

0. 

G 

294.15 

F 


refs react element 


Hydrngen(l) 

H 2. 


100. 

-2196. 

t 

20.2 7 

f 

.0709 

M (series D) 

14 (series D) 

Hydrogen peroxide 

H 2. 

C 2. 

100. 

-44*40. 

l 

294.15 

0 

1.407 

14 (series O 

14 (series C) 

IRFNA 

H \, 

*7216* 1.620650 4.69505F .C2699 

100. 

-64460. 

L 

294.15 

0 

1.44 

(•) 

<•) 

JP-5, AST MAI 

C t. 

M l.«||5 

100. 

-5103. 

L 

294. 1> 

F 

.407 

<0 

(0 

JP-4, RP-I 

C l. 

► 1.942J 

too. 

-44 JO. 

l 

298.15 

F 

. 771 

(R> 

X0 

Lithiunt(l) 

III. 


100. 

Wi4.1 

L 

959.8* 

F 

• 512 

• 

>0 

Lithium(s) 

III. 


100. 

0. 

5 

298.15 

F 

• 5 >4 

reference element 

so 

Lithium perchlorate 

Ill . 

CU. C 4. 

too. 

-90840. 

$ 

298.15 

0 

2.49 

14 (series A) 

14 (series A) 

Methanefe) 

C 1. 

t* 4. 

100. 

- 1 7894. 

G 

294.15 

f 


1 


Methane! 1) 

c t. 

* 4. 

IOO. 

- 2 1 990. 

l 

111. 44 

F 

.4219 

14 (series E) 

14 (series E) 

Methyl alcohol 

C 1. 

► 4. 01, 

100. 

-57040. 


298.14 

F 

.74499 

10 

31 

Monomethyl hydrazine 

C 1. 

► 4. 42. 

no. 

i 29C0 . 

l 

298.15 

F 

.474 

10 

SS 

Nitric acid 

H 1 . 

h|. 0 J. 

100. 

-41460. 

L 

298.15 

0 

1.9327 

14 (series C) 

14 (series O 

Nitrugcn(g) 

* 2. 


100. 

0.0 

C 

298.15 

F 


refs react element 


Nitr'ven(l) 

N 2. 


too. 

-2999. 

l 

77, >5 

F 

.404 

14 (series D) 

1 4 (series D) 

Nitrogen tetroxide 

a 2. 

C 4. 

100. 

- 46*0. 

L 

298.14 

0 

1.491 

14 (series Q 

14 (series Q 

Nitrogen trifiuoridc 

« 1. 

f 9. 

no. 

-J410J. 

L 

144.14 

0 

1.911 

14 (series E) 

14 (series E) 

Nitrome thane 

C 1. 

► 1* nu c 2 . 

too. 

-2 70 JO. 


294. | j 

F 

1.1171 

14 (series A) 

SO 

Octane 

C I* 

H IS. 

100. 

-5¥?40. 

l 

294.15 

f 

.49849 

XX 

XX 

Oitygen(g) 

C 2. 


1C0.0 

0.0 

C 

294.15 

0 


reference element 


Oaygen(l) 

0 2* 


no. 

- J l 02. 

t 

90.14 

n 

1.149 

1 4 (series D) 

14 (series D) 

Oxygen difluoride 

0 l. 

F 2. 

ICO. 

1069. 

l 

127.4* 

0 

1.921 

I, 14 (series E) 

14 (series E) 

Oxune{g) 

C 1. 


100. 

>4100. 

c 

294.15 

0 


0 


O tone!!) 

C 9. 


too. 

>0)10. 

l 

142.44 

c 

1.449 

14 (series E) 

14 (series E) 

Fentaborane 

• 5. 

h 4. 

100. 

7740. 

L 

294.1* 

f 

.4141 

14 (series A) 

14 (Miles A) 

Perchloryl fluoride 

cu* 

Cl. 7*. 

ICO. 

-11 >50. 

L 

224.4b 

0 

1.192 

10, 10, XI 

so 

Propane 

C ). 

M I. 

100. 

->0972. 

l 

211.0b 

F 

• 9404 

», sx 

XX 

n- propyl nitrate 

C 1. 

47. n i. ci. 

100. 

-51270. 


294.15 

9 

1.0994 

XI 

so 

Toluene 

C 7* 

4 S. 

100. 

2147. 

l 

294.15 

* 

.44210 

XX 

ss 

Unsym metrical dimethylhydratlne 

C 2. 

MS. N 2. 

100. 

*11900. 

l 

294.14 

F 

.741 

14 (series A) 

sx 


t, Mild; L, liquid; G, fU. 
b Fu*l, F; oaident, O. 

^Inhibited red fuming nitric add baaad on following weight p er c ents: HNOy(l) - M l, HjO^O) - 14, HjO(l) » X, HF(g) -0.1 
^Typical |tt fuel having following properties: H/C weight ratio - I. Ill , heat of combu s tion • II 000 Btu/lb. 

•Typical Jet tael having following properties: H/C weight ratio • 0. ID, heat of combustion • II 040 Bte/lb. 
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TABLE VIII. - VARIABLES IN INPT2 NAMELIST 


Variable 

Dimen- 

sion 

Type 

Common label 

Value 

before 

read 

Definition and comments 

KASE 

1 

, 

INDX 

0 

Optional assigned number associated with case 

P 

26 

R 

POINTS 

0 

Assigned pressures; chamber pressures for 
rocket problems; values in atm unless 
PSIA, NSQM or MMHG = T (see below) 

NSQM 

1 

L 


False 

2 

a Values in P array are in N/m 

PSIA 

1 

L 


False 

a Values in P array are in psia units 

MMHG 

1 

L 


False 

a Values in P array are in mm Hg units 

V 

26 

R 

POINTS 

0 

Volume, cm /g 

RHO 

26 

R 

POINTS b (P) 

0 

o 

Density, g/cm 

T 

26 

R 

POrNTS 

0 

Assigned temperature , K 

MIX 

15 

R 

MISC b (OXF) 

0 

Values of equivalence ratios if ERATIO = T; 
oxidant to fuel weight ratio if OF * T; per- 
cent fuel by weight if FPCT = T; and fuel to 
air weight ratio if FA = T 

eratio 

1 

L 

MISC 

False 

Equivalence ratios are given in Mix.® 

OF 

1 

L 

MISC 

False 

Oxidant to fuel weight ratios are given in MIX® 

FPCT 

1 

L 

MISC 

False 

Percent fuel by weight are given in MIX® 

FA 

1 

L 


False 

Fuel to air weight ratios are given in MIX® 

TRACE 

1 

R 

MISC 

0 

(5.E-9 for 
SHOCK 
problem) 

Option to print mole fractions 2 TRACE in 
special E -format 

IONS 

1 

L 

INDX 

False 

Consider ionic species® 

IDEBUG 

1 

I 

INDX 

0 

Print intermediate output for all points in- 
dexed 2: integer value 

TP 

1 

L 

INDX 

False 

Assigned temperature and pressure problem 

HP 

1 

L 

INDX 

False 

Assigned enthalpy and pressure problem® 

SP 

1 

L 

INDX 

False 

Assigned entropy (SO) and pressure problem® 

SO 

1 

R 

MISC 

0 

Assigned entropy, cal/(g)(K) 

TV 

1 

L 

INDX 

False 

Assigned temperature and volume (or density) 
problem® 

uv 

1 

L 

INDX 

False 

Assigned internal energy and volume (or den- 
sity) problem® 

sv 

1 

L 

INDX 

False 

Assigned entropy (SO) and volume (or density) 
problem® 

RKT 

1 

L 


False 

Rocket problem® 

DETN 

1 

L 


False 

Detonation problem® 

SHOCK 

1 

L 

INDX 

False 

Shock problem® 

SIUNIT 

1 

L 


False 

If true, the output tables will be in SI units 


®lf variable is set to be true. 

b Equivalenced to variable given in parentheses. 
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TABLE DC. - VARIABLES IN RKJTINP NAMELIST 4 


Variable 

Dimen- 

sion 

Type 

Common 

label 

Value 

before 

read 

Definition and comments 

EQL 

1 

L 

PERF 

TRUE 

Calculate rocket performance assuming 
equilibrium composition during expansion** 

FROZ 

1 

L 

PERF 

TRUE 

Calculate rocket performance assuming 
frozen composition during expansion** 

PCP 

22 

R 

PERF 

0 

Ratio of chamber pressure to exit pressure; 
list should not include values for the cham- 
ber and throat; storage allows for 22 values 

SUBAR 

13 

R 

PERF 

0 

Subsonic area ratios 

SUPAR 

13 

R 

PERF 

0 

Supersonic area ratios 

TCEST 

1 

R 


3800 

Optional initial chamber temperature esti- 
mate K. May be necessary when condensed 
species have been inserted on INSERT cards 
and 3800 K is far outside the raige of the 
data for the inserted species 

NFZ 

1 

I 


1 

Option for freezing composition at throat 
(NFZ * 2) or a supersonic point (NFZ > 2) 
when FROZ * . true. An extra table is 
printed with equilibrium composition 
through point NFZ and frozen thereafter 
with the composition of point NFZ as- 
Mmod. If NFZ >2, only 13-NFZ addi- 
tional exit points are allowed. 


^Required for rocket problems only. 

^et variable false If these calculations are not desired. 


TABLE X. - VARIABLES IN SHJONP NAMELIST 4 


Variable 

Dimen- 

sion 

Type 

Value 

before 

read 

Definition and comments 

INCDEQ 

1 

L 

TRUE 

Calculate incident shock parameters assuming 
equilibrium compositions** 

INCDFZ 

1 

L 

TRUE 

Calculate incident shock parameters assumir^; 
frozen compositions** 

REFLEQ 

1 

L 

FALSE 

Calculate reflected shock parameters assum- 
ing equilibrium composition 0 

REFLFZ 

1 

L 

FALSE 

Calculate reflected shock parameters assum- 
ing composition frozen at incident composi- 
tlon c 

U1 

13 

R 

0 

Shock velocity in m/sec (not required if values 
of MACH1 are listed) 

MACH1 

13 

R 

0 

Ratio of shock velocity to the velocity of sound 
in the unshocked gas (not required if values 
of U1 are listed) 


a Required for shock problems only . 

^t variable false if these calculations are not desired. 
c If variable is set to be true. 
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TABLE XI. - SAMPLE OF INTERMEDIATE OUTPUT. ITERATION DETAILS 


ITERATION 1 HAIKU 


0.775C006-01 0.3750006-01 0.3250006-01 0.175000E-01 0.525000E-01 

0.375C00E-QI 0 • 2 22 50''' *00 0 .6 750006-0 l 0. 5OO000E-02 O.R25000E-01 

0. 3?5C0Oc-0 1 0. <*>75 0006-0 l 0.2500006*00 0. 6250006-01 0.1100006*00 

C.I75C00F-01 0.5000006-02 0.6250006-01 0.180000E+00 0.7000006-01 

0. 525COOF-0 1 0 .8250006-0 1 0.1 100006*00 0.7000006-01 0. 

0 • 7*8 1 20 c *0 0 0 . 96082*6 *00 0.2252556*01 -0 . *326 066 -01 0.1295486*01 


0. 74 8 1 206*00 
0. 9608246*00 
0.2252546*01 
-0. 432686E-01 
0. 1295486*01 
0.3200916*02 


-0 . 1 4926 1 E *0 1 
-0 .25 3549 6 *0 1 
-3.2999476 *01 
-0.3602106 *01 
-0.2821706*01 
-0.2203276 *02 


PI N H C f 

-0.11*5976*02 -0.8 30 7886*0 1 -0.6709 166*01 -0.1638 776*02 0.229937E*00 0.2705426*00 

T« 0. 3800 J003E *04 6 NN ■ 0. 085999996*00 6NNL «-0 . 2302585 1 6 *0 1 PP» 0. 10000000E *0 3 LN P/N* 0.69077 5 52 6*0 1 0 .22 699 1 466 + 00 



NJ 

LN NJ 

DEL LN NJ 

CIST 

0. 

0. 

C. 1067566-15 

c 

0.250000E-02 

-0. 5991*66*01 

-0.1920866*00 

CF 

0.2500006-02 

-0.5991*66*01 

0.3759306*01 

CP2 

0.2S0C006-0? 

-C. 5991*66*01 

0. *5 701 16*01 

CF 3 

0.2500006-02 

-0. 5991*66*01 

0.2266*06*01 

CP* 

0.250000t-0? 

-0.59914,66*01 

0.72*8106*00 

CM 

0.2500006-02 

-C. 59914,66 *01 

0.65859*6-01 

CM2 

0 . 2500006-0 2 

-0.5991*66*01 

-0.8727606*00 

CM3 

0.25C000E-02 

-0.5991*66*01 

0. *899686*00 

CM* 

0.2500006-02 

-0.5991 *6E*01 

-0.1**6176*01 

CN 

0.2500006-02 

-0. 599 l*6p *01 

0.29*5696*01 

CN2 

0.250000E-02 

-0.59 9 1*66 * 0 1 

-0**235526*01 

C2 

0.2500006-0? 

-0. 5991*6E*0l 

-0.9138*36*00 

C2F 2 

0.2500006-0? 

-0.5991*66*01 

-0.99361 16*00 

C2P* 

0.2500006-02 

-0. 599 1 *6 e *01 

-0.282*056+01 

C2H 

0.2500006-0? 

-0. 5991*66 *01 

0.2696096*01 

C2MF 

0.2500006-02 

-0.5991*66*01 

0.200*516+01 

C2H2 

0.2500006-02 

-0.5991*66*01 

0.30*6066*01 

C2M* 

0.2500006-02 

-0,5991*66*01 

-0.276*576*01 

C2N 

0.2500006-02 

-0.5991*66*01 

0.20*2016+01 

C2N2 

0.2500006-02 

-0.5991*66*01 

0. 126559E+01 

C 3 

0 . 250000E -0 2 

-0.5991*66*01 

-0.1316106*01 

C* 

0.2500006-02 

-0.5991*66*01 

-0.6263116*01 

C5 

0.2500006-02 

-0.5991*66*01 

-0.71352*6+01 

F 

0.2500006-02 

-0.5991*66*01 

0.5075176*01 

FCN 

0.2500006-02 

-0.5991*66*01 

0.3032376*01 

P? 

0.250C006-0? 

-0.5991*66*01 

-0. 855*326 *00 

H 

0.2500006-0? 

-0.5991*66*01 

0.**3929E*0l 

HCN 

0.2500006-02 

-0. 5991*66*01 

0. *872*96*01 

MP 

0.2500006-02 

-0.5991*66*01 

0.8810676+01 

M2 

0.2500006-02 

-0. 5991*66*01 

0.5*16316*01 

N 

0 .250000E-02 

-0. 5991*66*01 

0.140*636-01 

NF 

0.2500006-0? 

-0.5991*66*01 

-0.4878576*00 

NF2 

0.2500006-02 

-0.5991*65*01 

-C. *039166*01 

Nf) 

0.2500006-0? 

-0. 599 l*6E *0 1 

-0.1030716*0? 

NM 

0.2500006-02 

-0.5991*66*01 

0.61*6206+00 

NH2 

0.2500006-02 

-0.5991*66*01 

0.2293296*00 

NM3 

0.2500006-02 

-0.5991*66*01 

-0.5664116*00 

N2 

0.2500006-02 

-C. 5991*66*01 

0.6631436*01 

N2C 

0.2500006-02 

-0. 5991*66*01 

-0.5190636+00 

N2M* 

0.2500006-02 

-0.5991*66+01 

-0.1008276*02 


HOJ/RT SOJ/R -COJ/RT -CJ/R7 

0. 25 56 706 * 0 1 0.6820846*01 0.4272 146*01 0 .4272146 *01 

0.2497871. *02 0 • 254243E *02 0 . 44 565 2F * 00 - C .4 7 C639E *0 / 

0.1215796*02 0.363812E*02 0.24223 , 6*02 0.2330706*02 

0.4242526*00 0.4499056*02 0.4*56636*02 C. 4365006*02 

-0.6191286*01 0.5421886*02 0.6041016*02 0.594s38t*02 

-0. 1827676 *02 0.6021916*0? 0. 78 495 HE ♦'•? C.775795F*02 

0. 2296806* 02 0.3252356*02 0.9555**6*^1 C. 8638146*01 

0.1821 37t*02 0.3642466*02 0. 1 8 21 l 06 *u2 0.17.9476*02 

0. 122663 1 *02 0.4175696*02 0.2949066*02 0.2857436*02 

0.7260306*01 0.4447706*02 0.3721676*02 C.363O0*E*Q2 

0.1794976*02 0.349345E*02 0.1698476*02 C.l6Ct85F*02 

0. 24B664t ♦ 02 0.4429846*02 0. 1 9432 06 *c 2 0.1851576*02 

0.3084146*02 0.3568836*0? 0.484694E*01 C.393065E*01 

0.9788096*01 0.5296656*0? 0.431784E*02 C.4?26?1E*02 

-0. 7005B5E*01 0.7160106*0? 0.7860686*0? 0.7769066*02 

0.2146936*02 0.4077046*02 0.193011E*C2 0.183846E*02 

0.1270506*02 0.5004266*02 0.373375E*02 0.3642126*02 

0.1583356*02 0.4531646*0? 0.2948296*02 C. 285*666*02 

0.1361356*0? 0.5464806*0? 0.408345E*02 0.3991826*02 

0.242775L*02 0.4535586*02 0.2107836*02 C .2^1 6206 *02 

0.18 72556*0? 0.5202916*02 0.3330366*02 C.3. !fc73E*02 

0.3143946*02 0.4243146*02 0.1059216*02 0.1007586*02 

0. 396707C *02 0.5019806*02 0.1052736*0? C. 9*1 1006*01 

0.4244256*02 0.5805696*02 0.1561446*02 C. 1 *69816 *02 

0.4037166*01 0.2564776*02 0.2081066*02 C . 1969436*02 

0. 759720C+01 0.4382716*02 0.362299E+02 C. 353 1366*02 

0. 4255241*01 0.3565046*0? 0.31 39516+02 0.30*7686 *02 

0.9203066*01 0.2014686*02 0.1094376*62 f.l0CV7*E*02 

0.1046056*02 0.3970616*02 0.2924566*02 0.2632936*02 

-0.4906166+01 0.3058386*02 0.3548996*02 C. 3*57376*02 

0. 3769146*01 0.2546796*02 0.21 69876+C2 C . 20 7 P2*e *02 

0.1727596*02 0.2480226*02 0.7526236*01 C .0 6C99* F *0 1 

0.1200676*02 0.3681436*02 0.2480766*02 0.2389136*02 

0.7542656*01 0.4636446*02 0.388217E+0? C . 37 5C5* 6 *02 

0.4659376*01 0.5435096*02 0.4969156*02 C. 4877526*0* 

0.1456896*02 0.3173626*02 0.171673E*02 C. 1625106*02 

0.1076336*02 0.3688256*02 0.2611926*02 0.252-296*02 

0.6120526*01 0.4100806*0? 0.348874E+0? C.3? C 711E*02 

0. 387959t*01 0.3314726*0? 0.292676E*(2 C. 28351 36*02 

0.2023006*02 0.4463206*02 0.24402C6+0? C. 235*576*02 

0.1531156*02 0.5800366*02 0.42*9226*0? 0.417759^*02 
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T 06RIV MATRIX 

0.2569696-01 

0.6660436-03 

0.190 388E-02 

0.6660436-03 

0.3841546-01 

0.1 89 197E-02 

0. 1983886-0? 

0.1891976-02 

0.696977E-02 

0.8654966-04 

0.3207386-01 

0.1753*16-02 

0.1367006-01 

0. 3684806-0 1 

0. *75*026-02 

PI N 

M 

C 

0.2271096+01 -0.4120886+00 

0.105*536+02 

P 06RIV MATRIX 

0.256969E-01 

0.6660436-03 

0.1903086-02 

0.6660436-03 

0.3841546-01 

0.1 89 197E-0 2 

0. 1983886-0? 

0.189197E-02 

0.6969776-02 

0.8654966-04 

0.3207386-01 

0.1753*16-02 

0. 1 367006-01 

0.368*806-01 

0 .4754026-0? 

PI K 

H 

C 

0.51 5C06E *00 

0.4*1 7*0E*00 

0.2799856*00 

POINT. 1 P. 

0. IOOOOOE*03 

T- 0.453*576 


M*« 0.2065816*02 CP/R- 0.6361396*00 


0.8654966-04 0.1367006-01 0.7157426*01 

0.3207386-01 0.3684806-01 -0.6786826-01 

0.1753416-02 0.4754026-02 0.7174556-01 

0.3826656-01 0. 3759726-01 -0.7883086-01 

0.3759726-01 0. 0.285781E-02 

F 

-0.1679276*01 -0.5330336*00 


0.8654966-04 0.1367006-01 0.1367006-01 
0.3207386-01 0.3684B0E-01 0.3684806-01 
0.1753416-02 0.4754026-02 0.4754026-02 
0.3626656-01 0.375972E-01 0.3759726-01 
0.375972E-01 0. 0. *860726-01 


F 

0. 631 931E *00 -0.3426846-01 
♦0* M/R- 0.1295906*02 S7R» 0.136009E*0 
OLVPT — 0. 1034276*01 OLVTP- 0. 1533036*01 


0. 

0. 

0. 

0. 

0. 


0. 

0 . 

0. 

0. 

0 . 


G AMH A ( S ) * 0.11 6900E ♦ C 1 V» C. 1801 1 76*03 
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TABLE Xn. - INTERMEDIATE OUTPUT, TEST FOR CONDENSED PHASES 


ALIS) 

300.000 

933.000 

IUSE" 

-l 

0. 




GO-SUMIAI J*Pt I) 

1 - 0.71701316*02 


MAX 

NEG 

DELTA 

G 

■ 0 . 

4111) 

933.000 

5000.000 

iuse* 

-1 

0 . 




ALCL3 ( S) 

300.000 

965.700 

IUSE" 

-2 

0 . 




ALCL3 ( L ) 

935.700 

3000.000 

IUSE" 

-2 

0 . 




G0-$UM( A I J9PI 1 » 

• 0. 3 785 7596 *02 


MAX 

MEG 

OEwTA 

G 

• 0. 

ALN(S) 

300.0)3 

3000.000 

IUSE* 

-3 

0 . 




GO-SUMUU*PI 1 ) 

* 3. 3 137253c +02 


MAX 

MEG 

OELTA 

G 

« 0. 

AL203 (SI 

300.000 

2315.000 

IUSE- 

9 

0.1667S09E-02 

AL203IU 2313.000 

5000.000 

IUSE* 

-9 

0 . 




CC SI 

300.000 

5000.000 

IUSE* 

-5 

0 . 




GO-SUMI A 1 J*PI I > 

• -0.13170536*01 


MAX 

MEG 

DELTA 

G 

• 0. 

H20 ( S ) 

200.000 

273.150 

IUSE" 

-6 

0 . 




H20CL) 

273.153 

373.150 

IUSE" 

-6 

0. 




MSISI 

300.000 

922.000 

IUSE* 

-7 

0 . 




GO-SUM(AU*Pin 

* 0.91375666*02 


MAX 

MEG 

DELTA 

G 

" -0.13170536*01 

M3<U 

922.000 

5000.900 

IUSE" 

-7 

0 . 




MGCL2 ( $) 

300.000 

987.000 

IUSE* 

-6 

0 . 




GO-SUM ( AI J*PI ( ) 

■ -0.10335976*02 


MAX 

MEG 

DELTA 

G 

■ -0.131 705 3 E * 0 L 

MGCL2ILI 

387.000 

>000.003 

IUSE" 

-8 

0 . 




MGOISl 

300.000 

3098.000 

IUSE" 

-9 

0 . 




GO-SUMI Al J*PI I ) 

- -0.13689916*02 


MAX 

NEG 

DELTA 

G 

■ -0. 10335WE*02 

MGOIU 3098.000 

5000.000 

IUSE- 

-9 

0 . 




S ( SI 

300.000 

388.360 

IUSE" 

-10 

0. 





S(U 388.1*0 2000.003 

G0-SUH< AI J*PI I ) - 0. 1 2911196*02 


IUSE* -10 0. 

H»< NEC OElT* G • -0.13689916*02 
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Figure L * Flow diagram for main routine. 





MAIis 


-REACT 

-SEARCH 

-THERM PH 


-SHCK- 


“NEVVOF — 

- SAVE 

- EQLBRM- 
“0UT1 

-0UT2- 

-OUT3 

-EQIBRM 

-OUT1 


-OUT3 

-CPHS 

-HCALC- 

-NEWOf- 

-SAVE 


-HCALC — 
-CPHS 
■MATRIX 
-GAUSS 

■EFMT 
- VARFMT 


CPHS 

MATRIX 

GAUSS 

“EFMT 

-VARFMT 


-CPHS 


-CPHS 

-HCALC- 


| EQLBRM 1 

OUT1 i 


-DETON- 


-OUT1 

-OUT2 

-OUT3 

-HCALC- 

“NEVVOF- 

“SAVE 


CPHS 

MATRIX 

GAUSS 

-EFMT 
■ VARFMT 


-CPHS 
-HCALC - 


i 

| EQLBRM— 1 “ 1 1 


-ROCKET- 


i 

E 


CPHS 
MATRIX 
GAUSS 

OUT1 

C 

OUT3 

VARFMT 

FROZEN CPHS 

SAVE 

NEWOF HCALC — 


-CPHS 


-CPHS 


-EFMT 

-VARFMT 


—CPHS 


Fioure Z - Subroutine tree oiaoram. 
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Figure 3. * Flow diagram of an application module. 
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